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00 (54) Title: SYSTEM AND METHOD FOR THREE-DIMENSIONAL IMAGE RENDERING AND ANALYSIS 

(57) Abstract: The present invention relates to methods and systems for conducting three-dimensional image analysis and diagnosis 
and possible treatment relating thereto. The invention includes methods of handling signals containing information (data) relating 
to three-dimensional representation of objects scanned by a scanning medium. The invention also includes methods of making and 
analyzing volumetric measurements and changes in volumetric measurements which can be used for the purpose of diagnosis and 
treatment 



o 



wo 01/78005 



PCT/USOl/11820 



SYSTEM AND METHOD FOR THREE-DIMENSIONAL IM ACT. RENDERING 

AND ANALYSIS 

BACKGROUND OF TH E INVENTION 

The invention was made with government support under R01CA63393 and 
5 ROl CA78905 by the National Cancer Institute. The government has certain rights to the 
invention. 

The present invention relates to the art of diagnostic imaging, and, in particular, to 
methods and systems for conducting highly enhanced three-dimensional reproduction and 
analysis of three-dimensional images aiid diagnosis and treatment related thereto, 

10 The present application claims priority to two (2) provisional applications 

identified as follows: U-S. Application Serial Number 60/196,208, filed April 1 1, 2000; 
and U.S. Application Serial Number 60/253,974, filed November 29, 2000. Each of 
these earlier filed provisional applications are incorporated herein by reference. 

It is known in the art of image diagnostic systems and methods to scan objects 

1 5 foxmd in a host patient, retrieve signals resulting fi-om such scanning, and process the 

signals to obtain information for possible diagnosis and treatment To date, much of the 
effort has been directed to treatment of such data in two dimensions, and analysis has to a 
certain extent been constrained by two-dimensional limitations. Consequently, 
professionals relying on such information have, similarly, been somewhat constrained in 

20 their ability to create diagnostic and possible treatment models. 

As a result of the inventors herein having the perspicacity to think and analyze 
based on three-dimensional vision, new methods and corresponding systems are now 
available for use in various applications. In one of the more basic applications, detection 
and analysis of objects in the body, especially growths, can be produced and analyzed to 

25 a degree of accuracy previously thought unknown. This is particularly useflil especially 
m the area of life threatening pathologies, e.g., cancer. One particularly pernicious 
pathology is lung cancer. The ability to be able to detect early stage lung cancer has 
provided strong motivation to improve diagnostic capabiUties. 
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For example, the small pulmonary nodule is a common radiographic finding that 
poses a diagnostic and clinical challenge in contemporary medicine. Approximately 
170,000 new lung cancers are detected each year in the United States alone. Pulmonary 
nodules are the major indicator of these malignant lung cancers, but may also be signs of 
S a variety of benign conditions. Of the all the new iimg cancers currently detected each 
year in the United States, approximately 90% are fatal, responsible for nearly 1 60,000 
deaths in 1999 alone [46], In fact, limg cancer is the number one cause of cancer 
mortality, responsible for more deaths in men than prostate and colon cancer combined, 
and in women, for approximately as many as breast and colon cancer combined [46]. 

10 One contributing factor for this low survival rate is the historical lack of early 

detection^ or the detection of the disease at a stage where it is most likely to be curable. 
Until now, the majority of lung cancers are either "symptom-detected," or found on chest 
x-rays (CXR) of asymptomatic patients. In most cases, the tumor is detected at an 
advanced stage where opportunity for successfvd intervention is limited. It has now been 

15 shown, through the work of the Comell Early Lung Cancer Action Project (ELCAP) that 
it may be possible to detect a significant number of malignancies at an early stage, where 
cure rates are much more promising, through the use of Ixmg cancer screening [28]. The 
comerstone of the ELCAP study is the use of computed tomography (CT) for the 
detection of small pulmonary nodules and the xxse of early-repeat CT (ERCT) for 

20 evaluation of these nodules, by high-resolution computed tomography (HRCT). It has 
been reported that low-dose CT has the potential to detect tumors 4-6 times as frequently 
as CXR [28]. It is theorized that this may engender a "stage-shift" in tumors detected by 
CT versus CXR. This additional time for potential intervention holds the promise of 
increasing survival from the current 12% to over 70% 5-year survival following resection 

25 of a stage / cancer [28, 76]. 

The tools developed in this work concern the application of computer vision and 
image analytic techniques to provide quantitative information describing the 
characteristics of pulmonary nodules seen in helical CT studies. They are an effort to 
capture in a reproducible, mathematical form the somewhat subjective characteristics 

30 traditionally used by radiologists and other physicians to describe these lesions in order to 

2 
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determine their likelihood of malignancy. Advances in helical CT technology are now 
allowing for the detection of a significant number of pulmonary nodules less than 1 cm in 
diameter. At this point, traditional methods for nodule characterization (many based on 
hand-measurement and subjective evaluation) become less effective, given the small size 
5 of these lesions. Through the use of computer-aided diagnostic techniques, we have a 
unique opportunity to assist the radiologist in the increasingly difficult task of diagnosing 
these small nodules. Furthermore, as the expected widespread acceptance of CT 
screening for early detection of lung cancer increases, computer-aided tools should play a 
prominent role in the early detection, characterization, and cure of lung cancer in many 

1 0 thousands of patients around the world. 

Pulmonary nodules, while most important when they represent malignancy, are 
firequently the manifestations of a variety of benign conditions such as hamartoma, 
benign adenoma, or rarely, sarcoidosis, amyloidosis or infection. Fimgal pulmonary 
infections (histoplasmosis, cryptococcosis, coccidiomycosis), mycobacterial 

15 (tuberculosis) and other bacterial infections may also, infrequently mimic neoplasms 
[69]. 

Establishing a specl&c diagnosis is important, since surgery is usually not needed 
in the treatment of benign small puhnonary nodules (SPNs). On the other hand, SPNs 
may represent malignancy (cancer) of several types (squamous-cell, adenocarcinoma, 

20 large-cell, small-cell, bronchioloalveolar carcinoma, etc.) where expeditious removal by 
thoracotomy is usually indicated- For example, approximately 5% of patients with 
"small-cell" limg cancer are diagnosed after a chest radiograph (x-ray) shows an SPN. 
Currently, the definitive diagnosis is made after the tumor has been removed surgically. 
If the diagnosis is made before the operation, clinical decision making is facilitated since 

25 mxxltiorgan scanning and inediastinoscopy to detect spread or metastasis are performed 
prior to thoracotomy [85]. In addition, very small puhnonary nodules may not be 
detectable in the ordinary chest x-ray or conventional CT and are visible in heUcal 
(spiral) CT examinations. Therefore, improvement in the available techniques for the 
detection (e.g. differentiation from vessels) and differentiating malignant from benign 

30 SPNs are needed and will result in clinical and economic benefits. 
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Small pulmonary nodules are usually identified on routine chest radiographs 
(CXR) in asymptomatic patients or on chest radiographs obtained for other indications. 
Most nodules have a benign cause, but bronchogenic carcinoma, metastases, and other 
malignant processes are common causes and must be excluded in the differential 
5 diagnosis. As they are potentially malignant, SPNs require expeditious evaluation. 

Techniques including computerized tomography (CT) have been used for the differential 
diagnosis of pulmonary nodules. Radiologic evaluation includes chest x-ray (CXR) or 
computed tomography (CT) for the detection of nodules, high-resolution CT (HRCT) for 
further characterization, and CT and magnetic resonance imaging (MRI) to evaluate the 

10 mediastinum, and extrathoracic structures. 

The radiologic characteristics that differentiate benign firom malignant pulmonary 
lesions include size, the presence of calcification in a "benign pattern," the shape and 
surface characteristics, rate of change in size^ clinical characteristics of the patient such as 
age and smoking history, etc. Meticulous assessment of the margin of the nodule is 

1 5 necessary. Useful characteristics for the differentiation of benign JBrom malignant 
include, sharpness of margin, the presence of spicules, length of spicules, spicules 
extending to the visceral pleura, pleural tail sign and circumscribed pleural thickening. 
Even with current criteria however, the best sensitivity and specificity that can be 
achieved are 85% and 78%, respectively. Further improvements in diagnostic accuracy 

20 are needed since at the present time many patients with ambiguous findings are referred 
for early surgery [72]. 

Mediastinal CT is the preferred modality for examining the mediastiniun in these 
patients while magnetic resonance imaging is used selectively, e.g. in patients Avith 
superior sulcus tumors who are candidates for surgery. The differential diagnosis is an 

25 important consideration, particularly in patients with small pulmonary nodules, due to the 
increased sensitivity in lesion detection, increased specificity and lesion (tissue) 
characterization (usually done with MILI) [29]. Metabolic imaging by positron-emission 
tomography (PET) provides many advantages over conventional anatomically-based 
cross-sectional imaging, although its use is generally resolution-limited to nodules larger 

30 than 1 cm in diameter. For routine use in oncology, a detailed assessment of specific 

4 
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efficiency of PET is indicated [67], Mediastinoscopy is done using CT for guidance [61]. 
Fine-needle aspiration biopsy is a useful diagnostic tool [95]. Thoracic biopsy is a safe 
and accurate minimally invasive surgical approach to resectable peripheral SPN. Many 
patients require surgical biopsy especially those with adequate physiologic reserve (that 
5 is those who would be able to tolerate surgery) [53]. Also, video-assisted thoracic surgery 
(V ATS) often allows one to accomplish the same goal as the comparable open procedure 
but with less morbidity and a shorter hospital stay. With continued development of 
instrumentation, increasingly complex procedures continue to be accompUshed with this 
technique [37]. 

1 0 Currently the diagnosis and management of SPNs is based on the basic principle 

that every nodule must be regarded as potentially malignant until proven otherwise. 
Malignant nodules are surgically removed unless the patient is in such bad shape from 
other diseases that he or she can not survive the surgery or there is proof that the cancer 
(malignant tumor) originated elsewhere in the body and that the SPN is a metastasis. On 

1 5 the other hand, resection of a benign nodule carries a small surgical cost The major 
benefit derived from resection of the benign SPNs is that the diagnosis is made by 
pathologic^ examination and malignancy is excluded. In addition, once a small 
pulmonary nodxile has been detected, the diagnosis should be made quickly to avoid 
spreading of the malignancy. 

20 When a malignant cause cannot be ruled out, the patient's age, smoking history, 

and nodule size must be considered. For moderate and high-risk patients, an immediate 
and more invasive work-up is indicated [1 7]. Observation by serial radiographs or CT 
studies may be the appropriate course for patients who are at low risk for malignancy. 
When the probability of malignancy is low, the patient may be advi§ed to have repeated 

25 radiologic examinations m order to determine whether the SPN is changing over time 
before the decision about whether to operate is made. Therefore, making the diagnosis 
without surgery either at initial examination or after consideration of repeated 
examinations would be of benefit [47]. 

In addition to primary puhnonary pathology, pulmonary nodules may be due to 

30 metastatic disease from cancers of the breast, stomach, pancreas, kidney, bladder, and the 

5 
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genitoxarinary tract. Computed tomography, especially helical CT, is probably the most 
sensitive imaging technique in the identification of pulmonary metastases, because it 
detects a higher number of pulmonary nodules compared with other techniques. In a 
population of patients presenting Avith multiple, solitary pulmonary nodules, or the 
5 absence of nodiiles, helical CT shows 30 to 40% of supplementary nodules when 

compared to conventional scaiming techniques. HeHcal, or spiral, scanning is a relatively 
new computed tomographic technique based on the continuous acquisition of volumetric 
CT data during continuous x-ray beam rotation and continuoxxs patient translation through 
the scanner at constant velocity. It has many advantages over conventional CT including 

10 rapid image acquisition during a single-breath-hold scan of the lung, and the ability to 
obtain axial image reconstruct ions at arbitrary and overlapping intervals, thus allowing 
the detection of small lesions that otherwise would be inconspicuous because of 
respiratory motion or slice averaging. This leads to better identification of small 
pulmonary nodules and to high quaUty multiplanar reconstructions that can be useful in 

1 5 the study of mediastinal lymph nodes and the vascular and tracheobronchial spread of 
lung cancer [10]. 

In addition to its use in standard radiographic studies, helical CT also enables the 
optimal use of contrast products in which the iodine is poorly concentrated to study the 
pulmonary vessels. The continuity, both anatomically and for the lesions obtained by 

20 spiral CT is such that it is now possible to apply to thoracic pathology techniques of 

multiplanar and three-dimensional reconstruction. If all the information is contained in 
conventional transverse imaging is shown slice by sUce, not everything is perceived by 
the observer because the information is inconvenientiy presented, by deconstructing 
anatomical and lesional picture. In helical CT, reconstructions of tiie volume may be 

25 inspected from a coronal, sagittal, oblique or three-dimensional perspective, furnishing 
supplementary information to axial images alone [66]. 

Traditionally, measurements of pulmonary nodules have been made by the 
radiologist using calipers or a ruler on a chest radiograph. The goal was to make an 
approximate measure of the nodxile size by estimating the diameter seen in the projection 

30 captured on the chest film. The sizes could be recorded, along with subsequent measures 
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to establish whether the nodule was growing and at what rate. This practice has carried 
over into the reading of n[iany CT examinations. Although computer-based workstations 
are sometimes available for the review of CT studies, the sequential images are often 
output m "hard copy," film format and viewed on a Ught box. Here too, caliper 
measurements are the norm, with size being estimated from the image with the greatest 
cross-sectional area. 

As a step toward precision and reproducibility of measurement, many digital 
review workstations provide "digital calipers," that allow the radiologist to make 
measurements using a mouse. While a general improvement of caliper-based 
measurements on film, it should be noted that the estimate of nodxde size is still observer- 
dependent, for the following reasons: (a) selection of an appropriate cross-section 
(maximal area) for measurement is variable; (b) nodule cross-sections may not be 
circular, prompting multiple measurements and averaging schemes; (c) the margin of 
many nodules is somewhat indeterminate; and (d) window and level settings can 
markedly change the apparent size of a nodule. 

These factors suggest that a more automated assessment of nodule size would be 
beneficial to the unbiased, reproducible estimation of size and growth. A potentially 
more important point is that the measurements made by the radiologist are nearly always 
two-dimensional in nature. With computer-aided methods, the entire nodule volume may 
be measured with great precision, eliminating the problems of slice selection, non- 
circularity, and window and level distortion. 

Computer vision is the study of techniques for the extraction of information 
fi:om images. The images may be two-, three-, or higher-dimensional in nature, and are 
acquired by varying types of image sensors, including conventional 2D cameras, video 
microscopy systems, x-ray film digitization, and computed tomography systems. The 
underlying physics of acquisition also varies, including the absorption of visible light, 
infrared, x-ray, and the sensing of magnetic resonance. Each of these may be 
accomplished in a variety of geometries and projections. 

Computer vision is closely alhed with its sister disciplines: image processing 
and computer graphics. Lnage processing is the appUcation of techniques for the 

7 
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improvement or modification of an image, generally so that it may be better perceived 
by a himian observer. Computer graphics is concemed with basically the inverse 
problem of computer vision: generation of images based on an informational 
description. 

5 Most of the techniques now popular in computer vision have been developed for 

industrial and military applications. Examples include product quality control, and 
detection of tactical targets in radar, sonar, and conventional images. More recentiy, 
computer vision has been applied to problems in biology and medicine. Automated 
analysis of cell morphology from microscopic images, analysis of anatomical structures 

10 radiographic images, and study of function using advanced imaging modalities are fast 
becoming key applications for computer vision methods. 

One notable difference between computer vision in medicine and biology and 
traditional computer vision is that in industrial applications, the structure of objects being 
studied is frequentiy known a priori, and also often has a carefully described geometry. 

15 Biological structures, however, while they normally have a particular structure, exhibit a 
wide range of variation, especially in pathological cases. Even normal anatomical 
structures defy simple geometric description, as they are three-dimensional, complex 
structures that exhibit variation between subjects and also often have a djmamic 
appearance in the same subject. Growth over time, as well as motion during image 

20 acquisition make it non-trivial to obtain easily-modeled data of a particular structure, 
even in a single subject. 

These subtleties motivate techniques that can characterize the size, shape, and 
other qualities (e.g. density) of a biological structure in a quantitative, reproducible way. 
This allows for the comparison of data between subjects, as well as the study of a 

25 structure in a single subject over time. Such measures are precisely what is needed to 
effectively study the progression of a pulmonary nodule from first detection to eventual 
diagnosis, and in many cases, d\iring follow-up as well. 

The purpose of computed tomography (CT) is to display the anatomy of a slice of 
the body by measuring absorption of x-rays as they pass through the body in many 

30 trajectories. Imaging usually is done in slices perpendicular to the axial (head to toe) 

8 
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direction. The resulting measurements are manipulated mathematically to obtain two- or 
three-dimensional displays. 

X-rays are electromagnetic radiation (10'^^ - 10'^ nm) that are produced when an 
electron beam strikes a metal anode, typically made of tungsten. The anode is heated 
5 during the generation of the x-rays and dissipation of heat becomes a limitation in long 
examinations. Only x-rays in one general direction are allowed to exit the x-ray 
generator by the use of a beam collimator made of lead shielding material. As they travel 
through matter, the x-rays are attenuated, i.e., their intensity (related to the number of 
photons per irnit cross sectional area) is decreased by absorption or scatter. Absorption 

10 occurs through the photoelectric effect, where the x-ray energy is given to an electron 
tiiat is freed from the atom. Scatter occurs where x-rays of higher energy (generated by 
electrons accelerated by higher voltages, having higher frequency) transfer some of their 
energy to an electron (the energy necessary to free it from the atom), while the remaining 
energy becomes a secondary photon (x-ray) that is transmitted in a new direction. Thus, 

1 5 absorption is useful in diagnostic radiology since it occurs in a known direction while 
scatter creates noise due to secondary x-rays of unknown direction being detected by the 
x-ray sensor(s). 

CT and all other x-ray diagnostic examinations are based on the fact that different 
structures of the body absorb or scatter x-rays to a different degree. The frmdamental 
20 principle is expressed by Beer's Law 

§=-^ (1.1) 

aL 

describing the attenuation of intensity of electromagnetic radiation as it passes through a 
homogeneous medium. In this equation, / is the intensity of the radiation (number of x- 
rays, photons, per unit surface area); L is the length of the pathway, and yu, is the linear 
25 attenuation coefficient. This equation merely describes that the x-ray attenuation as it 
goes through a imit length of the medium is proportional to the intensity of the incident 
radiation or that each x-ray (photon) has equal probability of being absorbed while 
traveling through a given thickness of a given material as any other x-ray (photon). This 

9 
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differential equation leads to the following exponential relation describing the intensity / 
of radiation as it passes through a medium: 

I^he'^ (1.2) 

where / is the transnutted intensity through a thickness L, This is derived easily from the 
5 fundamental definition of e 

4-(.e') = e' (1.3) 
ax 

and the associated chain-rule for differentiable functions of x 



-{e-he-- (1.4) 



10 



Therefore, 



dL dL ^ov y ^i.j; 

15 Figure 1 illustrates Equation 1 .2, where the intensity of the radiation passing 

through the material, /, is a function of the intensity of the incident radiation, /o> the 
length of thematerial, and the linear attenuation coefficient of the material, //. This 
coefficient is determined by the density of the material, the atomic number, and the 
wavelength (or energy spectrum for polychromatic x-rays) of the x-ray radiation imder 

20 consideration. Thus, knowing the linear attenuation coefficient in each voxel of the body 
would provide information on the density and atomic mmiber of the elements at that 
location. 

The process of creating images from x-ray absorption data along a variety of 
angles and translations is called image reconstmction, involving the solution of a set of 
25 equations relating the geometry of absorption data seen in each projection. A given CT 
projection is the result of x-ray traversal through objects of varying attenuation 
coefficients at varying distances. This relationship can be generalized as 

10 
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h = V - MiA • V^^^^ -he - /x„i„ (1.6) 

where / (/„) is the intensity of the radiation emerging jfirom the body and detected by the 
detector, fiu 112,^^- iWj, are attenuation coefficients of the individual structures traversed by 
the x-ray and Zy, £2. . - • ^5. are corresponding distances traveled by the x-rays (lengths). 
Assuming equal lengths firom x-ray source to detectors in helical CT (L, chosen at the 
time of image reconstruction), Equation 1 .6 can be rewritten as 

K ^I,e-^^'^^'''-^^''^L (1.7) 

Figure 2 illustrates the physical system in which Equation 1 .7 applies. In this 
case, the incident x-ray beam passes through several materials of varying length 
and attenuation characteristics. The computation of the individual attenuation 
coefficients (fiiji2^...,pi^ requires the solution of an equal number of equations obtained 
fi'om the different projections acquired in the CT procedure. 

The usual method for solution of the set of equations and the generation of the CT 
image is cB\\Qd filtered back-projection, based on the Radon transform. In this method, a 
uniform value of attenuation is projected over the path of each ray such that the total 
attenuation is proportional to the measured attenuation [14, 15, 6], These values are 
stored as elements of one vector of the reconstructed matrix. The procedure is repeated 
for each ray sxmi of the CT while corrections are made for voxels traversed by the x-ray 
beam in an oblique fashion. The assumption that x-ray attenuation is uniform throughout 
the path of each x-ray beam (obviously an oversimpUfication) results in blurring of the 
resulting image. In spite of this limitation, filtered back-propagation is popular because of 
its inherent parallelism. It allows processing of the data obtained by each ray while data 
acquisition continues in other projections, dramatically improving computational 
efficiency. The blurring is decreased by increasing the number of projections and by 
convolution with a filter function (convolution kemel). The filtering functions depend on 
x-ray tube geometry, detectors, intended effects (sharpening edges thus enhancing spatial 
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resolution at the expense of decreased density resolution versus more refined density 
resolution while sacrificing spatial sharpness, etc.). 

The relative pixel values obtained by this image reconstruction process are not the 
absolute values of the linear attenuation coefficients of individual voxels of the body (/^i 
5 +M2 Mn)' The CT numbers, called Hounsfeld Units (HU) in honor of Godfirey 

Hoxinsfield (the inventor of CT scaiming [30, 31]), are related to the linear atteniiation 
coefficients as follows: 

CT number = ^^^"^"^ (1.8) 

where /ia, is the attenuation coefficient of water, ju the attenuation coefficient of the pixel 
10 in question, and K a constant. The value of K must be large enough (e.g. 200 - 2000) to 
accommodate the accuracy of the scanner. For example, consider a CT scanner with a 
density resolution of ±0.5%, or ±1 in 200. In this case, a value ofK= 200 would be 
sufficientiy high to encode density values, as they are typically recorded with integer 
precision. A larger value of would expand the scale beyond the accuracy of the 
1 5 scanner. 

Helical CT is a newer method of data acquisition in computerized tomography 
[83, 40]. Data acquisition is continuous during rotation of the x-ray source and the 
detector (mounted in a toroidal gantry) around the patient who is simultaneously moved 
in the axial (head to toe) direction through the rotating gantry at precise speeds from 1 
20 mm to more than 10 mm per second. The resulting helical projections are used to 

reconstruct an anisotropic Cartesian space representing the x-ray absorption at each voxel 
in the scanned patient volume. 

In commonly used instruments, the x-ray beam is coUimated to a fan that defines 
the image plane at any given time, and the array of detectors travels in a circular path 
25 aro\md (typically 360°) the patient. This approach provides the opportunity to obtain 
many projections (and therefore images) in a short period of time (e.g. in one breath- 
hold), and the collection of data in a continuous volumetric manner, rather than in slices. 
This allows reconstruction in any position as well as improved two-dimensional and 

12 
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three-dimensional views. The major practical advantage of spiral CT is the ability to 
cover complete anatomic regions in a single exposure [38]. Isotropic imaging, where 
spatial resolution is the same in all three directions, will eventually become a reality [39]. 
Helical CT has been found useful in the study of pulmonary nodules where optimally 
centered studies have been shown to improve accuracy [40, 88]. 

The most important parameters of a CT protocol relate to the image geometry, 
reconstruction resolution, and radiation dose. Pitch is defined as the ratio of the table 
speed in mm/revolution (assuming a 360'' reconstmction algorithm) to the collimation 
width of the x-ray beam. For a given beam collimation, the table may be advanced at a 
speed that would move the table through one collimation width per revolution, or 1 : 1 
pitch, or at a higher pitch, reducing the effective radiation dose and scan time at the cost 
of image quality. For example, a protocol with 10 mm collimation and 20 mm/revolution 
table speed would acquire images with a 2: 1 pitch. 

The reconstruction resolution is described by several parameters. The in-plane 
resolution (the image plane is perpendicular to the scanner axis), which is equal in the x 
and y dimensions, gives the resolution of a pixel in a single 2D reconstructed image. The 
axial resolution, or resolution in the z dimension, is determined by the slice spacing, or 
the physical distance between reconstructed 2D images. This spacing is a function of the 
scan collimation, pitch, and reconstruction algorithm. Furthermore, images may be 
reconstructed at overlapping intervals (e.g. 10 mm overlapping slices at 5 mm intervals). 
The axial resolution of conventional CT protocols is often 5-20 times more coarse than 
the in-plane resolution. This results in 3D voxels that are anisotropic (not all dimensions 
exhibit equal spacing). 

Radiation dose is a function of x-ray tube current, typically expressed in mA, tube 
voltage, expressed in kV (peak measiirements are given units of kVp), and scan duration. 
The tube current and voltage directly affect the quality of reconstructed images, as the 
signal-to-noise ratio seen at the x-ray detectors is proportional to the dose used. 
Improvements in detector sensitivity and filtered image reconstruction algorithms are 
continually reducing the minimum dose required to achieve acceptable image quality at a 
given scan geometry. 
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There are three types of CT examinations commonly performed in the study of 
puhnonary nodxiles. The first is a low-dose, low (axial) resolution screening scan, aimed 
at the detection of possible pulmonary nodules in the lungs. The second, following 
detection of a pulmonary nodule, is a standard-dose scan of the entire lung volume. Hie 
5 third (and primary focus of this work) is a standard-dose, high-resolution scan of the 
region containing the detected pulmonary nodule, used for analysis and characterization 
of the nodule itself. A description of the parameters typically xised in the different scan 
protocols is shown in Figure 3. 

Current CT scaimers impose a tradeoff between image quahty, sUce thickness, 

1 0 and x-ray dose. Many scanners have a gantry rotation speed of approximately 1 

revolution/second. This implies that a breath hold on the order of 30 seconds to 1 minute 
is required to obtain a whole lung scan. More recent scanners, using multislice 
technology, have reduced this time to just a few seconds, helping to eliminate respiratory 
and other motion artifacts. It may be anticipated that future developments in CT detector 

1 5 : and processing technology may make possible much more rapid and higher resolution 
scans such that future screening scans may be taken at high («0-l rmn) axial resolution. 
However, current CAD systems are designed for the constraints of today's scanners 
which suggest a two-resolution approach: a low-dose initial whole-lung scan followed (in 
some protocols) by a high-resolution focused scan of detected nodules. 

20 Much research has been done relating to the manual detection and 

characterization of pulmonary nodules by radiologists. In the past 10-15 years, work in 
computer-assisted and automated detection has appeared in the Uterature. Until quite 
recently, the dominant imaging modality for the study of cancer in the lung has been the 
chest radiograph, or chest x-ray (CXR). With the increasing use of CT for the detection 

25 of these lesions, however, we are now seeing an increase in the use of automated methods 
for the detection of pulmonary nodules in thoracic CT studies. Some work has also been 
reported dealing with characterization of pulmonary nodules, but it is only quite recently 
that advances have been made in this area. As characterization of pulmonary nodules 
firom CT scans is the focus of this work, the following sections provide a view of 
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pertinent related issues in the study of nodules using both CXR, and CT, in the screening 
setting, and for nodule characterization. 

The goal of lung cancer screening is to detect pulmonary nodules, establish their 
probability of malignancy at an early stage, and thereby improve likelihood of cure. The 
S traditional detection of pulmonary nodules occurs by their appearance on a routine chest 
radiograph, or one taken for unrelated reasons. In the past 10-15 years, research has been 
done in developing computer-aided methods for detection of these nodules. 

In the 1970s, the National Cancer Institute (NCI) sponsored three large screening 
trials to evaluate the benefits of lung cancer screening using chest radiography and 

10 sputum cytology. These were held at the Mayo Clinic, Johns Hopkins, and Memorial 
Sloan-Kettering. [52, 20, 19]. Although these studies found that lung cancer screening 
allowed for earlier detection of lung cancer, better resectability of tumors, and better 
survivorship of the surgery, the overall mortaUty in the screened and control groups were 
similar. Thus, lung cancer screening was not recommended as national policy. 

1 5 Detection of puhnonary nodules has traditionally been done through visual 

inspection of chest radiographs. Since the 1970's, however, much research has been 
devoted to the automation of this process. Computer-aided detection and analysis 
systems for pulmonary nodules fall into several categories. Image processing techniques 
have been used to improve contrast and otherwise enhance candidate regions for the 

20 visual detection of nodules [43, 44, 77]. Systems have been developed for the automated 
detection of candidate nodules [50, 96, 93, 8, 59]. Some work has also been done on 
automated analysis of nodules in CXR images for the prediction of malignancy [71, 58, 
87]. 

The preprocessing stage of several CAD systems for the detection of nodules in 
25 CXR involves histogram manipulation. Techniques have been explored using histogram 
equalization and high-frequency enhancement [59]. Several studies have described the 
use of image differencing, which is the subtraction of a "nodule suppressed" image from 
one that is "nodule enhanced," to help reduce the influence of other anatomical structures 
appearing in the radiograph [23, 96, 44, 93]. This enhancement of CXR images is 
30 achieved using one or more matched filters designed to identify the round regions 
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characteristic of many nodules. For example, a disk-shaped filter may be convolved with 
the image to enhance nodule-like regions. A ring-shaped median filter may be used for 
nodule suppression. Related work in preliminary segmentation of lung images has also 
been done in the areas of lung boundary delineation [91,3] and rib identification [70]. 
5 These steps help eliminate regions of the image fi'om consideration as potential nodules. 

A second stage of histogram manipulation may be performed following the initial 
image preprocessing. Multiple thresholding of the preliminary images may be performed 
to identify the incremental contribution of each intensity band to candidate regions. In 
this way, suspect regions are iteratively enlarged (by decreasing a threshold) and 

10 analyzed [50, 8]. Computed tomography (CT) is fast becoming the modality of choice 
for the detection and analysis of pulmonary nodules. Studies in lung cancer screening 
with CT have been performed in Japan [41, 78] and New York [28]. 

There have been several prototype systems for performing complete, automated^- 
Ivmg CAD [94, 22, 84, 7] (i.e., taking as input whole-limg scans and identifyiug nodules). 

15 Note that these systems involve low-resolution screening scans and therefore, may be 
able to detect small nodules (<1 cm) but there is insufficient information to characterize 
them. These early attempts employ a multiphase approach involving the following three 
stages: (a) identification of the lung regions in the CT images; (b)separation of candidate 
nodules fi"om vessels within the lung regions; and (c) classification of candidate nodules. 

20 For diagnostic classification of small nodules, a focused HRCT scan of the 

region-of-interest is obtained. Prototype CAD systems for these scans have been 
developed in which the region-of-interest is specified by the radiologist and just the 
separation and classification stages are automatic [63, 42]. 

The lungs may be separated firom other structures in the chest by first performing 

25 some noise reduction through spatial filtering and then applying a simple threshold with 
backgroimd removal to the whole-lung scan. As an example, a single slice from a whole- 
lung screening scan is shown in Figure 4 and the extracted lirng region is shown in Figure 
5. The entire lung volume can be measured or visualized as shown in Figiore 6. A similar 
method was used by Giger et al. [22], while Brown et al. [71] used an anatomical model 

30 with fuzzy set membership based on several soft constraints including threshold, 
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continuity and relative position. Toshioka et al. [84] use geometric filtering to avoid 
irregularities in the detected boundary. They also use a second method, basing their lung 
boundary on predictions from rib locations in order to identify the presence of very large 
nodules. 

5 Further analysis of other thoracic structures is also possible. For example, the 

lumen of the trachea and the main bronchi may be traced by a simple seeded region 
growing algorithm; typical results for screening data is shown in Figure 7. This figure 
illustrates the liroitations of the approach on screening scans due to their low axial 
resolution; the bronchial tree and the vascular structure can be extracted to a much finer 

10 degree with higher-resolution scans. Three-dimensional region growing combined with 
tree-traversal algorithms has also been used in removal of the lung vasculature as a 
preliminary step to nodule detection [90, 16, 79]. 

Until recently, nearly all of the literature on the characterization of pulmonary 
nodules from radiologic images was concerned with the radiographic appearance of the 

1 S lesions to a trained radiologist. Most of the measures of nodule size were made manually 
using calipers or a ruler on chest radiographs or hardcopy CT images. Measurements of 
shape and density distribution were also described in somewhat qualitative terms. 

Some notable work on the characterization and subsequent classification of 
pulmonary nodules may be found in [26, 47, 82, 73]. The chief characteristics reported 

20 were size, smoothness of edge, presence and (for larger nodules) pattem of classification, 
spiculation, lobidation, cavitation, and the presence of a pleural tail. 

Perhaps more important than shape and density characterization, the estimation of 
nodide aggressiveness as a fimction of growth rate has been studied for some forty years. 
Nodide doubling time, as computed from manual measures of nodule diameter, has long 

25 been utilized effective predictor of nodule malignancy [12, 56], 

A variety of techniques have been used to characterize potential candidate regions 
for nodule detection in screening data Measurement of object circularity is often helpfril, 
as this metric may help exclude rib crossings and other structures [50, 44, 92]. Methods 
involving gradient analysis have also been described [50, 49]. These are based on the 
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notion that the edge gradients of a small nodule-like object should primarily point toward 
the center of that object, unlike the case with ribs and many vessels. 

Morphological operations have been explored for the detection of nodules in CXR 
[96]. These nonlinear filters impose geometric constraints on the candidate regions, 
5 allowing, for example, round objects (e.g. nodules) to be distinguished from longer, 
flatter ones (e.g. ribs, heart, etc.). Correlation-based template matching has also been 
used to identify candidate nodules, an idea related to matched filter enhancement [8, 59]. 

The selection of candidate regions may be refined by one or more automated 
decision-making steps. Both neural networks [93, 59] and rule-based [8] systems have 
10 been used in this capacity. Inputs to these decision-refinement stages come from one or 
more shape, edge, or histogram features. 

One of the most important issues regarding the use of CAD systems for nodule 
detection is the effective detection rate. Most computer-aided schemes for the analysis of 
chest radiographs unfortunately have relatively low sensitivity. In addition there is 
1 5 firequentiy a high false-positive rate (very low specificity) in these systems. A typical 
result is 2-7 false-positives per image, with sujB&cient sensitivity (over 70%) [93, 8]. 
Thus, the choice of operating point on the receiver operating characteristic (ROC) curve 
for a particular CXR analysis technique may be quite challenging. 

In addition to the detection of piilmonary nodules, computer techniques have 
20 recently been studied for the classification of nodules found in chest radiographs. 

Computation of density gradients as a measure of shape and surface irregtdarity has been 
used to distinguish benign from malignant nodules [71]. Fractal texture analysis has also 
been applied to the nodule classification problem with encouraging results [58, 87]. 

When assessing nodule growth for the prediction of malignancy, two or more CT 
25 studies, separated by a suitable InterScan interval, are needed. It is theorized that it may 
also be possible to estimate the probabiUty of nodule malignancy based on size, shape, 
and density parameters determined firom a single exam. 

In a screening setting, 2D features are firequently measured to refine the set of 
nodule candidates for detection. In this capacity, Giger et al. [22] used measures 
30 including perimeter, area, compactness, circiilarity, elongation, and location within the 
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lung. Similarly, Toshioka et al. [84] used area, thickness, circxilarity, density mean and 
variation, and location. For small nodule characterization, however, HRCT studies are 
required (see Figure 3) and methods for high-resolution characterization are being 
explored [42, 63, 51]. 

5 McNitt-Gray et al. [5 1] have explored two-dimensional (2D) sh^e metrics 

including aspect ratios, circularity, and compactness, as well as a variety of texture 
measures (e.g. correlation texture and difference entropy). Their study also employed 
bending energy as a measure of surface characteristics. Three-dimensional measures of 
surface characteristics were studied by Kawata et al. (42). Using data isotropically 

10 resampled to 0.33 mm^ voxels, they measured Gaussian and mean curvature as metrics of 
irregularity in nodule margin. 

Another challenging problem in the study of pulmonary nodules is the analysis of 
non-solid nodules, or ground-glass opacities (GGOs) [18, 1 1]. These lesions are difficult 
to characterize using traditional size and shape metrics. 

IS Patient demographic information has been used as an adjunct to radiographic 

analysis in the prediction of small pulmonary nodule status. Several patient 
characteristics have been shown to be highly correlated with nodule malignancy. These 
include advanced age and smoking history (pack-years) [81, 82, 47]. 

To study the pulmonary nodule using computer-aided techniques, we must 

20 formulate one or more models of what a pulmonary nodule is, in particular, with respect 
to its radiographic appearance. Pulmonary nodules are lung neoplasms that fall into two 
general categories: benign and malignant. The separation between these categories is 
generally one based on rate of growth. Most benign lesions do not grow, or grow much 
more slowly than malignancies. There are, however, some notable exceptions. Infections 

25 in the lung may exhibit extremely fast growth. Also, it has been suggested that there may 
be some indolent, slow-growing malignant species. These malignancies are by far the 
exception, however, and doubling times for most maUgnant nodules range firom 30-400 
days [47]. 

Pulmonary nodules appear radio graphically as round, opaque lesions with density 
30 slightly more dense than that of water (-- 0-1 00 HU). Their position in the lung and 
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immediate surroimdings, from an image-analytic perspective, differentiates them into one 
of the following categories: (a) well-circumscribed - the nodule is located centrally in the 
limg, without significant connections to vasculature; (b) vascularized - the nodule is 
located centrally in the lung, but has significant vascularization (connections to 
5 neighboring vessels); (c) pleural tail - the nodule is near the pleural surface, connected by 
a thin structure ("pleural tail"); and (d) juxtapleural - a significant amoiint of the nodule 
periphery is connected to the pleural surface. 

Hie above categories describe the types of solid pulmonary nodules seen in CT 
images. A second class, known currently as ground glass opacities (GGO) have a 
1 0 significantly less uniform radiographic appearance and will be considered separately later 
in this work. 

One of the best predictors of nodule malignancy is growth rate. Similarly, nodule 
size, is also highly correlated widi malignancy [82]. Significantly larger nodules (> 3 cm 
in diameter) are more likely represent lung cancer [82]. Overall size has therefore been 
1 5 used in characterization and classification studies [26, 82]. 

In addition to these benefits, volumetric measurement allows for the identification 
of anisotropic growth. It is commonly assumed that puhnonary nodules grow uniformly 
in all dimensions. 

In 1956, Collins et al. [12] succinctly described the relationship of growth and 
20 cancer, as well as the need for quantitative meeisurements of tumor growth: 

^*The definition of cancer, its diagnosis and its prognosis all depend upon description of 
growth To the layman a synonym for cancer is a 'growth ' There are no quantitative 
terms for the description of growth or growth rate in clinical use, " 
Their seminal paper continues, providing some of the most important early discussion of 
25 tumor "doubling time," a measure of the amount of time (usually expressed in days) for 
a tumor to double in volimie. Nodtde growth, and the rate at which it takes place, is 
perhaps the feature most predictive of malignancy. Nodule doubling time and 
calcification in "benign" patterns have been the only universally accepted diagnostic 
tools of thoracic radiologists in the differentiation of benign from malignant lesions. 
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Traditional measures of the size and shape of puhnonary nodules were made 
manually on chest radiographs using calipers. Nodule size was characterized by .a 
measure of diampter (based on the assumption that nodules appeared roughly circular), 
or as a combination of several **diameter" measurements. 
5 Assessment of pulmonary nodule shape has been based on subjective shape 

characteristics includin g notions of roundness, spiculation, lobulation^ and sharpness of 
edge. Nodule surface characteristics are frequently used to help classify benign from 
malignant nodules. In particular, the presence of spiculation on the surface, as well as 
lobulation of the overall nodule shape, have been reported as more common in 

10 malignant nodules [47, 82]. As with measurement of nodule size, there have been 

difficulties in making reproducible shape characterizations, due to the lack of a precise 
mathematical basis and measurement technique. 

Several techniques for measuring nodule surface characteristics have been 
described in the literature. Huo et al. [32, 33] used analysis of radial edge gradients (a 

1 5 technique originally developed for mammography) in two-dimensional CT images. 

Kawata et al. [42] used a method of curvature estimation based on 3D image gradient 
analysis. Such techniques for curvature estimation, and the prerequisite 3D gradient 
analysis and edge detection have been described in detail [5, 54]. 

Detection of pulmonary nodules in screening CT scans is challenging due to the 

20 desire to limit radiation dose and scan time, both of which result in a reduction of 

information. The need to scan the full lung volume in a single breath hold motivates 
the use of large intervals between images (10 mm). An additional reduction in 
radiation dose achieved through low tube current further affects image quality. With 
these constraints, computer algorithms for the detection of pulmonary nodules have 

25 traditionally focused on two-dimensional (2D) algorithms or serial 2D methods that 
examine coimectivity between nodule candidates in adjacent slices. 

Automated CAD systems for tixe detection of pulmonary nodules have been 
based on a multistage approach, as illustrated in Figure 8. The main stages are: (a) 
preprocessing, (b) identification, and (c) classification. Preprocessing involves noise 

30 filtering and isolation of the lung region from the image data. Identification is the use 
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of image filters and knowledge based reasoning to develop a set of nodule candidates 
from the lung volume being searched. Classification is the analysis of nodule 
candidates, normally as part of a rule-based decision system, to reject other structures 
(e.g. vessels) from consideration and to assess the likelihood that each remaining 
5 candidate is a true nodule. The classification stage in detection systems is sometimes 
extended to include the separation of benign from malignant lesions, however, this 
may only achieve limited results with the low-dose, low (axial) resolution studies 
common to initial screening exams. Characterization and classification techniques that 
lead to prediction of mahgnancy are generally perfomied using high-resolution CT 
10 (HRCT). 

A number of automated CAD systems for the detection of pulmonary nodules in 
low-dose screening studies have been described in the literature [94, 22, 84, 7, 4]. 

A typical lung cancer screening protocol consists of a low-dose (140 kVp, 40 mA) 
CT scan of the full lung volume using 1 0 mm coUimation, reconstructed at 5 mm 

1 5 intervals, with an in-plane resolution of 0.5 mm or better. Following reconstruction, the 
resultant voxels are highly anisotropic (0.5 x 0.5 x 10 mm). Although the in-plane 
resolution may be sufficient for the detection (and preliminary characterization) of many 
lesions, small pulmonary nodules (less than 1 0 mm in diameter) may appear in only one 
to three contiguous images. This anisotropy presents special challenges in distinguishing 

20 between small pulmonary nodxiles and fine vasculature. 

Blood vessels may have either a linear or elUptical cross-section in a single CT 
image, depending on orientation. The more perpendicular to the image plane a vessel 
lies, the more circular the resultant image. Therefore, it is a non-trivial problem to 
distinguish between small pxilmonary nodules and small vascular cross-sections as the 

25 low axial resolution of screening studies limits the three-dimensional information needed 
to fully characterize these objects. The problem is further exacerbated by the fact that 
nodules are frequentiy attached to the pulmonary vasculature or to the pleural surface, 
increasing their difficulty of detection as they may be considered part of either 
confounding structure due to the similar density values found in many non-calcified 

30 nodules. 
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The preprocessing stage identifies the region of the CT images that correspond to 
the lung volvune. Algorithms have been developed to automatically detect the thoracic 
cavity, lung volume, tracheobronchial tree, and pulmonary vasculature. 

The first step in most pulmonary CAD systems is to identify the lung boundary, 
5 thereby eliminating other thoracic structures firom the nodule search region. This is 
typically accomplished using serial 2D gray-level thresholding techniques on each CT 
image. Spatial filtering may be applied to reduce the effect of noise on the boundary 
segmentation. Once a preliminary boundary has been determined, geometric filtering 
may also be used to eliminate discontinuities in the boimdaries determined in adjacent 

10 sHces, as vs^ell as to reduce the contribution of image artifacts. 

An example of this procediure is shown in Figures 4 through Figure 6. Figure 4 
shows a single 2D image from a low-dose screening scan. Figure 5 shows the lung 
region in this cross-section produced using the techniques described above. A three- 
dimensional surface-shaded rendering can also be produced, as is shown in Figure 6. 

1 5 Giger et al. [22] used histogram analysis to determine appropriate thresholds for 

lung segmentation. Toshioka et al. [84] used spline representation and curvature analysis 
to refine lung boundaries based on standard thresholding as well as the detection of rib 
locations. This enables the detection of juxta-pleural nodules, which can be excluded 
firom the lung volume using simpler threshold-based methods. Anatomical models were 

20 developed by Brown et al. [7] as the basis for a knowledge-based system that utilized 

fuzzy set membership based on density, spatial locality, and position. These models were 
then used in the determination of lung, tracheobronchial tree, and thoracic cavity regions. 
In our studies, we have found that simple thresholding combine with geometric filtering 
has produced good results in identifying these regions. 

25 Another preprocessing step in pulmonary CAD systems is the identification of the 

tracheobronchial tree and pulmonary vasculature. Segmentation of the airways may be 
accomplished by seeded region-growing techniques. An example of this algorithm 
appUed to screening data is shown in Figure 12. It is evident that the low axial resolution 
of anisotropic screening data limits the extent to which the bronchi can be followed. 

30 Tree-based algorithms that explicitly track the bifiircatipns of the bronchi and pxilmonary 
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vasculature have been described by Wood et al. [90], Related methods based on 
morphological filtering [97] and fuzzy logic [103] as well as the benefit of airway 
subtraction to nodule detection have also been described [16]. 

Once preprocessing has eliminated regions outside the lung firom consideration, 
5 the next step in nodule detection is the identification of nodule candidates from the 

remaining lung volume. This is typically done through the use of one or more gray-level 
thresholding operations, as the nodules are of higher density than the surrounding lung 
parenchyma, Yamamoto et al, [94] combined single and double-thresholding with 
connected-component analysis to identify nodule candidates that are within a selected 

10 density range and size. In addition, morphological filters, including disk and ring-shaped 
kemels as well as 2D morphological opening ("rolling ball filter"), were used to select 
candidates. Multiple thresholding based on histogram analysis has been used to select 
nodule candidates and differentiate them from vessels [22] . In this system, while 
increasing, the threshold, vasculature that bifurcates along the scanner axis will produce 

1 5 one or more "daughter nodes" in adjacent slices, thereby reducing the likelihood that the 
candidate is a tme nodule. 

The final stage of a pulmonary nodule detection system, classification, is the 
refinement of a set of nodule candidates based on size, shape, and other measures. A 
likelihood estimate that the remaining candidates are nodules is also often assigned. The 

20 techniques used in refining the set of nodule candidates are somewhat related to those 
used in the characterization and classification of nodules in high-resolution diagnostic 
studies. The difference, however, is that the aim of nodule candidate classification is to 
separate nodules from non-nodules, while characterization and classification in diagnostic 
HRCT is used to separate benign from malignant nodules. With this distinction in mind, 

25 many of the nodule feature metrics used in detection and diagnostic CAD systems 

describe similar size and shape characteristics. It is important to note, however, that the 
standard low-dose screening study is significantly anisotropic, and therefore largely 
limits the feasibility of performing measurements in three dimensions. For this reason, 
features iised for nodule candidate classification have been primarily limited to 2D 

30 metrics or analysis of connectivity between adjacent 2D slices, 
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When a pulmonary nodule is detected at initial screening, high-resolution 
computed tomographic (HRCT) images may subsequently be acquired to further study 
the lesion. These studies typically involve scanning the region-of-interest (ROI) 
containing the nodule at a high in-plane resolution (0.3-O.S mm), with small beam 
5 collimation (1 .0 mm) and 1 : 1 pitch. At the time of HRCT, a diagnostic CT scan of the 
chest is also commonly done for use in patient management. Figure 3 shows a 
comparison between low-dose screening, diagnostic, and HRCT studies. 

An important consideration for HRCT studies is that the entire 3D volume of the 
lesion should be acquired without significant artifact (due to respiration or other patient 
10 motion). In traditional high-resolution studies, only the several images containing the 
largest 2D nodule cross-section were considered vital, as measurements would be made 
visually or with calipers. Thus, it is an object of the present invention to provide the 
capability for enhanced size, shape, and densitometric measurement of objects found in 
the body. 

15 SUMMARY OF THE INVENTION 



The present invention includes methods, systems, and computer readable 
medium for producing and analyzing three-dhnensional images. In particular, the 
present invention relates to image representation and analysis for diagnosis and possible 
20 treatment of humans. 

With respect to image production, the present invention includes providing a 
three-dimensional image representation of an object which detectably attenuates a 
signal resulting firom scanning. A signal corresponding to the three-dimensional image 
is obtained and can be supersampled in three dimensions to a degree greater than the 
25 highest spatial firequency of the signal. 

Scanning-medium used herein can include x-rays, magnetic field (as provided 
by MRI), ultrasound, laser, electrical current, visible light, iiltraviolet light, infrared 
light, and radio firequency. A preferred embodiment of the present invention 
contemplates the use of x-rays such that the signals produced therefirom provides CT 
30 tomography for processing and analysis. 
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Preferably, supersampling includes imposing three-dimensional units of image 
representation (voxels) on a three-dimensional image signal and increasing the number 
of such units in order to reduce the contribution of error per xmit. Preferably, the 
supersampling is conducted in three-dimensional isotropic space. 

This method can be conducted on a host body, especially a patient, wherein the 
object is a tissue mass such as a growth, e.g., a tumor. It has been foimd especially 
iisefiil for detecting and analyzing pulmonary nodules, and, especially, small pulmonary 
nodules which are eqxial to or less than 10 nun in diameter. 

After the image has been supersampled, it can be further processed by 
segmenting the object from other structures detected in the signal retrieved from the 
scan. Such further processing can include subtracting attachments such as vascular 
attachments from the object, and, subtracting the object from adjacent structures such as 
pleural structures. The image can be further processed by restoring volimietric 
characteristics of the object. 

Another aspect of the present invention includes three-dimensional 
segmentation of a region of interest and a host body to differentiate an object found in 
the region. The region of interest is scanned with, for example, x-ray, to provide a 
signal corresponding to a three-dimensional representation of the region upon which 
three-dimensional units of representation are imposed. The units are then identified as 
having object signal, backgroxmd signal and as boundary units which have both object 
and background signal. The image representation imits are then subjected to 
thresholdiug in order to separate background from object, followed by three- 
dimensional connected component analysis to identify those images as objects 
contiguous ia three-dimensional space. 

The connected component analysis can include component labeling in three 
dimensions such as by recursive connected component labeling or iterative cotmected 
component labeling. 

Furthermore, the segmentation procedure can include the step of three- 
dimensional supersampling of the three-dimensional signal image. 
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Segmenting can also include morphologically filtering images identified as a 
result of the steps set forth above. Such filtering can include three-dimensional opening 
of the images where by small connections between regions are broken and sharp edges 
are removed, such opening including »osion followed by dilation of the images in three 
5 dimensions. The method can also include a closing of the three-dimensional region of 
the background signal resulting from the previous steps, such closing including dilation 
followed by erosion of the background in three dimensions. 

The segmentation can also include three-dimensionally regrowing images to 
approximate true surface characteristic of the images. In one embodiment the 
10 regrowing is iterative constrained dilation in three dimensions. 

The segmentation procedure can also include eliminating structures adjacent 
images of interest, and, in the case of a nodule, such structures can include thoracic 
structures such as a pleural siurface. In order to remove thoracic structures, the present 
invention contemplates determming in three dimensions angles describing orientation 
15 of the surface of a stracture to be eliminated followed by performing an opening 

operation (based on the angles previously found) to detect a majority of the structure to 
the exclusion of the image of interest. Thereafter, the stmcture is three-dimensionally 
subtracted firom the image of interest. In a preferred embodiment, the step of 
determining in three dimensions the angles describing orientation of the surface of the 
20 structure can be conducted by three-dimensional moment analysis. Moreover, the 

method for eliminating adjacent structure can also include morphologically filtering in 
three dimensions the image of interest. Finally, the signal can be processed to provide a 
smooth surface representation of the image of interest. 

A smooth surface representation of a three-dimensional voxel representation of 
25 a three-dimensional image can include: segmenting the image to provide a segmented 
voxel three-dimensional signal; followed by modified tessellating the three-dimensional 
image to provide a three-dimensional triangulated curve surface estimation; after which 
the triangulated curve surface estimation can be smoothed by filtering. 

Another aspect of the present invention is the ability to provide an estimated 
30 volumetric change time (DT) of an object which imdergoes change in size in the body, 
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such as a growth or even an organ found in a mammalian host, e.g., a human. This 
method contemplates obtaining a measurement of a change of a volimietric 
characteristic of the object over a period of time. The "change" referred to herein can 
be growth or regression qf growth. The change of the volumetric characteristic is then 
5 related to the time period of such change by an exponential change model whereby an 
estimated doubling time can be derived by comparing the volxmie change over such 
time period. An "exponential change model" as used herein includes "exponential 
growth model." Thus, when the term "exponential growth model" is used in the 
description, exponential change model is deemed to be implicit. 

10 This method is implemented by obtaining two three-dimensional volxmietric 

assessments of an existing object at a first time ti and at a later second time t2^ The 
volumetric assessment can be a volume Vi followed by measurement of a volume F2. 
Alternatively, the volumetric characteristic can be a measure of the mean density at 
time ti and a measure of the mean density at time fe. Another embodiment 

15 contemplates combing a change retardation factor, iJ, with the standard exponential 
growth model. 

In a preferred method of estimating the doubling time, the mcEisurements can be 
obtained from three-dimensional images provided from signals having three- 
dimensional data resulting from scaiming the host, and, in other preferred embodiments, 

20 when such image has been supersampled, segmented, and both supersampled and 
segmented. This method is particularly useftil in conducting diagnosis and possible 
treatment of pulmonary nodules, especially small pulmonary nodules which are usually 
not greater than about 10 mm in diameter. 

An alternative embodiment of this method and system include the ability to 

25 estimate doubling time of a an object which does not appear in a first scan of the region 
of interest but does appear at a second scan at time /2* hi this case a minimum 
detectable change dimension is identified followed by calculating doubling time using 
a volumetric measurement at t2 and a volume measurement derived from the minimal 
detectable change dimension over the time period between ti and tz. 
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Another embodiment of the present mvention related to time between scans 
includes a method and system for determining time between scanning examinations, or 

interscan intervals. In this aspect of the invention^ a maximum magnitude of error, €, 

between volumetric measurements from scatming and actual object size is developed. 
5 A nninimuTTi reliably-detectable percent volume change (a) is then identified based on 
the maximum magnitude of error. A relationship is established between the magnitude 
of error and object size, preferably a linear relationship, and an expression using an 
exponential change model is derived in terms of change in volume and time. Based on 
the expression resulting from the exponential change model relationship, a time interval 

10 to observe a percent change is determined. Preferably, this method and system is 
highly applicable to potentially cancerous grovsrths, especially when the interscan 
interval is sufficiently short in duration to increase probability of early detection as 
cancerous. Once again this procedure is especially helpful in detection, analysis, and 
diagnosis of small pulmonary nodules, e.g., those not greater than 1 0 mm in diameter. 

1 5 The present invention also contemplates a method and system of determining 

volumetric change in size over a period of time. In this method a doubling time of 
change is first determined, followed by establishing a relationship based on a 
exponential change model of said period of time between volumetric change time and a 
relative change in volume, followed by determimng the volumetric change over such 

20 period of time. 

Another method and system of analyzing three-dimensional scanned-produced 
images includes three-dimensionally registering a region of interest taken at separate 
times. In one method the registration is effected by locating the center of mass (COM) 
in the region of interest and aligning subsequent images of the region of interest for 

25 analysis. This can be accomplished by iteratively determining the center of mass as 
explained in detail hereinafter. 

In another manifestation of this embodiment, the three-dimensional registration 
of regions of interest is conducted by three-dimensional correlation. In this case two 
images of region of interest are selected wherein the limits of search in such region are 
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specified The second region of interest is translated to all locations in the search 
region and the three-dimension correlation between the two regions of interest are 
determined to find the greater value of a match metric. 

Analysis of three-dimensional registration of two regions of interest can also 
5 include a process of weighting such as by use of Gaussian window, a triangular radial 
window, or a rectangular radial window, which can be used to reduce unwanted 
confounding structures found in the region of interest The three-dimensional 
weighting of the region of interest can also include method of moment analysis to 
provide density distribution of the weighted region of interest 

10 Yet another aspect of the present invention is a method of estimating curvature 

of an object which has been seamed to provide a signal having a three-dimensional 
representation of such object This method includes providing a three-dimensional 
triangularly tessellated representation of the object followed by determining the surface 
normal to all triangles in the tessellated representation. In this regard it is noted that 

1 5 every plane surface has a normal, thus each triangle has a surface normal which is 
determined in accordance with this method. Next, a surface normal at each vertex is 
then calculated based on the surface normals resulting from the previous step. Then the 
augular differences between the surface normal at each vertex (designated home vertex, 
Vi) and the vertex normals of all vertices adjacent to the home vertex are determined. 

20 Finally, the curvature is estimated at each vertex of the object based on the angular 
differences resulting from the previous step. The relationships and algorithms 
supporting this method have been set forth in detail in the Detailed Description of the 
invention. 

The present invention also includes a method of providing a cvirvature metric 
25 analysis of a three-dimensional object which has been scaimed to provide a single 
corresponding to such object. In this method, a measurable three-dimensional 
representation of an object is provided which includes a surface. Surface curvature is 
estimated of such surface at a mmiber of discreet points on the surface to provide a 
number of surface curvature estimates. Then a frequency distribution of the surface 
30 curvature estimates resulting from the previous step is determined. Subsequent to 
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determining the frequency distribution of the surface curvature estimates, various 
aspects of the distribution can be detennined, including, but not limited to variance, 
range of distribution, mean of the distribution, minimum of the distribution, and 
maximum of the distribution. Preferably the method set forth above can be applied to 
5 an object which is a growth in a human such as a pulmonary nodule, especially one 
which is not greater than 10 mm in diameter. 

Moreover, the analysis based on three-dimensional registration can also include 
estimated doubling time of an object based on a change in the mean density of the 
region of interest . In this case, an object such as an object image, is Gaussian- 
10 weighted and the estimation can be an iterative technique to select an appropriate 

standard deviation, a , for weighting. The doubling time has been calculated based on 

standardized mean density and a stopping criteria, e. 

As a result of the present invention, three-dimensional resampling techniques 
have been developed and evaluated which mitigate partial volume effects in the 

1 5 assessment of object volume. There has been a significant improvement in the 
reproducibility of volxmietric estimates afforded by the techniques of the present 
invention over use of standard reconstructed CT data. For example, the double time 
estimates are consistent Avith pathologic diagnosis and are more appropriate estimates 
than those based on 2D metrics. Furthermore, as a result of estimating curvature based 

20 on a smoothed tessellated polygonal surface model, the actual nodule surface is 
approximated given the available resolution. 

Another benefit of the present invention is to establish new techniques for 
computer-aided diagnosis (CAD) of small pulmonary nodules (SPNs). A collection of 
three-dimensional size, shape, and density metrics have been developed for 

25 characterization of small pulmonary nodules. This development has enhanced the 
ability to differentiate benign from malignant lesions. 

Furthermore, new methods have been developed for possible computer-^ded 
analysis and diagnosis of small pulmonary nodules in computed tomographic data. 
SmaU pulmonary nodxiles, those with a diameter of less than 1 0 mm in diameter, are 
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now commonly detected and analyzed according to their size, shape, density and 
growth characteristics. 

For a better understanding of the present invention, together with other and 
further objects, references made to the following description taken in conjxmction with 
the accompanying drawings and the scope of the invention will be pointed out in 
claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 is an illustration of Beer's Law in a homogeneous material; 
Figure 2 is an illustration of Beer*s Law in a heterogeneous material; 
Figure 3 illustrates a comparison of screening and diagnostic CT protocols; 
Figure 4 illustrates a single slice from a whole-lung CT scan; 
Figure 5 illustrates the lung region identified from the CT slice shown in Figure 4; 
Figure 6 illustrates a 3D reconstruction of the lung s\irface; 
Figure 7 illustrates a 3D reconstruction of the trachea and main bronchi; 
Figure 8 illustrates a multi-stage approach for the automated detection of 
pulmonary nodules in screening scans; 

Figure 9 illustrates a well-circumscribed small pulmunary nodule; 

Figure 10 illustrates a small pulmonary nodule with confoxmding vasculature; 

Figure 1 1 illustrates a small pulmonary nodule with a pleural tail. 

Figure 12 illustrates a small juxtapleural nodule; 

Figure 13 illustrates a model of a vascularized nodule; 

Figure 14 illustrates a volume overestimation at point of vessel attachment; 

Figure 1 5 illustrates a volume of a spherical cap; 

Figure 16 illustrates a volume overestimation of an 8 mm spherical nodule with 
vascular attachments of varying sizes following opening with a spherical kernel; 
Figure 17 illustrates a model of a jvixtapleural nodule; 

Figure 18 illustrates a two-dimensional sampling of a circle where n equals 0; 
Figure 19 illustrates a two-dimensional sampling of a segment of a circle where n 
equals 1 and 2; 
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Figure 20 illustrates a two-dimensional area estimation of a segment of a circle; 
Figure 21 illustrates a two-dimensional image of a small pulmonary nodule; 
Figure 22 illustrates a histogram of gray-level intensities in nodule image; 
Figure 23 illustrates a two-dimensional segmentation of a small pulmonary nodule 
5 using iterative threshold selection; 

Figure 24 illustrates segmentation results by varying threshold from —800 HU to 
100 HU; 

Figure 25 illustrates examples of morphological operations using a spherical 

kernel; 

10 Figure 26 illustrates a 2D morphological opening on a pulmonary nodule; 

Figure 27 illustrates a three-dimensional morphological opening for vascular 
subtraction; 

Figure 28 illustrates a three-dimensional iterative morphological filtering for 
vascular subtraction; 

15 Figure 29 illustrates a three-dimensional iterative morphological filtering for 

vascular subtraction. 

Figure 30 illustrates segmentation of a synthetic vascularized nodule; 
Figure 3 1 illustrates sensitivity of vascular subtraction as a function of vessel 
diameter and kernel size; 
20 Figure 32 illustrates sensitivity of vascular subtraction to kernel size in 

segmenting in vivo pulmonary nodules; 

Figure 33 illustrates an example of the pleural surface removal technique; 
Figure 34 illustrates an example of the pleural surface removal technique; 
Figure 35 illustrates a 15 conical 2x2x2 neighborhoods and associated 
25 triangulations in the marching cubes algorithm; 

Figure 36 illustrates a comparison of 3D surface representations of a known 

sphere; 

Figure 37 illustrates a comparison of exponential nodule growth for three 
doubling times TD = 75 (aggressively malignant), = 150 (malignant), and TD = 500 
30 (benign); 
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Figure 38 illustrates RMS percent error in 35 volume measurements of spherical 
silicone nodules; 

Figure 39 illustrates CT images of a small nodule at baseline and 33 days later; 

Figure 40 illustrates 3D visualization of nodxile at baseline and 33 days later; 

Figure 41 illustrates two-dimensional images and segmentation of a malignant 
small pxilmonary nodule at baseline and 36 days later; 

Figure 42 illustrates 3D visualization of malignant nodule exhibiting non-uniform 
growth; 

Figure 43 illustrates an ellipsoid with its principal axes labeled; 
Figure 44 illustrates a silhouette of an object whose shape will be discretized by 
the sampling procedure; 

Figure 45 illustrates the discretized representation of the object shown in Figure 

44; 

Figure 46 illustrates a polyline representation of the object shown in Figure 44; 

Figure 47 illustrates a filtered (smooth) polyline representation of the object 
shown m Figure 44; 

Figure 48 illustrates a comparison of 3D surface representations of a small 
pulmonary nodule; 

Figure 49 is an illustrative patch of a 3D tessellated surface; 

Figure 50 illustrates a computation of the surface normal of each triangle; 

Figure 51 illustrates a calculation of the surface normal at each vertex; 

Figure 52 illustrates surface normals forming the basis of curvature estimation at a 

vertex; 

Figure 53 illustrates an atomic measure of cvirvature based on the angular 
difference between vertex surface normals; 

Figure 54 illustrates a reduction in variance of curvature estimates of synthetic 
spheres when using average normal curvature over mean curvature; 

Figure 55 illiistrates a surface curvature analysis of a malignant small pulmonary 

nodule; 

Figure 56 illustrates a mean of average normal curvature distribution in synthetic 
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spheres; 

Figure 57 illustrates a spherical perturbation model for surface cxirvature analysis; 
Figure 58 illustrates a frequency distribution of curvature estimates in a synthetic 

sphere; 

5 Figure 59 illustrates a firequencey distribution of curvature estimates in surface 

perturbations; 

Figure 60 illustrates a frequency distribution of curvature estimates of a 
spherically-perturbed synthetic sphere; 

Figure 61 illustrates a frequency distribution of curvature estimates of a 
10 spherically-perturbed synthetic sphere; 

Figure 62 illustrates a variance of curvature estimates as a function of smoothing; 

Figure 63 illustrates a Receiver-Operating Characteristic (ROC) curve for the 
complete set; 

Figure 64 illustrates a Receiver-Operating Characteristic (ROC) curve using 
1 5 cro ss- validation; 

Figure 65. illustrates an iterative center of mass determination in a CT scan of a 
small pulmonary nodule; and 

Figure 66 illustrates a correlation-based registration of nodule ROIs. 

20 DETAILED DESCRIPTION OF THE INVENTION 

The present invention relates to methods and systems for conducting three- 
dimensional image analysis and diagnosis and possible treatment modalities relating 
thereto. For the purpose of disclosing the present invention and the present preferred 
embodiments thereof, the detailed description has been broken down by subject matter in 
25 order to more clearly explicate the various aspects of the invention. Headings have been 
provided for the sake of convenience and focus of subject matter, but are not meant in 
any way to limit the scope of applicability of the subject matter found in one section. 
Indeed, each technical aspect is generally considered germane to subject matter found in 
other sections. Thus, one must read all of the sections as related to the extent that the 

30 technology disclosed therein can be applied to the technology found in other sections. 
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Fuifhennore in this regard, a certoia degree of redundancy has been built into the 
detailed description in order to provide as much as possible overview and correlation of 
the technology between sections foxind herein, 

A system in accordance with one embodiment of the present invention can 
5 include a scanner, a processor, a memory, a display device, input devices, such as a 
mouse and keyboard, and an internal bus connecting the various components together. 
The system can be coupled to a commimication medium, such as a modem connected to 
a phone line or a wireless network, or an Internet connection. Since the components of 
a computer and how they operate are well known to one of ordinary skill in the art, they 

10 will not be discussed here. 

A scaimer includes any device capable of imposing a scaiming medium on 
subject, e.g., a patient, to produce a signal attenuated by objects in the field of the scan. 
These devices include, but are not limited to, CT scanners, MRI, PET, ultrasoimd, 
optical tomography, and electrical impedance tomography. A preferred scanner is a . 

1 5 helical CT scanner which was used in the work and analysis reported herein. 

A processor as referred to herein is a system of one or more electronic 
components which are statically or dynamically configured (in hardware or software) to 
implement an algorithm or process. Examples include systems based on 
microprocessors, application-specific integrated circuits (ASICs), field-programmable 

20 gate arrays (FPGAs), and/or discrete components. Any type of storage device, such as a 
RAM or ROM, can be used as the memory in the system. Additionally, the memory 
can be a magnetic, optical, magneto-optical, or solid state drive that is coupled to the 
processor and which can receive portable memory devices, such as floppy disks, hard 
disks, optical disks, or solid-state memory on which the method of processing computed 

25 tomography images in accordance with the present invention is programmed. All media 
which can store information, operating instructions, etc., and can be read by a computer 
for processing, storage, etc. is referred to collectively herein as "machine readable- 
medium." The method of processing computed tomography images in accordance with 
the present invention can be stored as an executable program in the memory. 
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Section 1: Computer Assisted Diagnosis (CAD) -Virtually all aspects of the 
present invention can be used to provide computer assistance in the process of diagnosis 
for detecting, treating, and monitoring the biological condition(s) of a patient. The 
degree to \sdiich the techniques and technology disclosed herein are used and/or relied on 
5 depends on the condition, the patient, the predictability of outcome, the physician and 
other factors too niunerous to detail herein. To the extent the technology disclosed herein 
can be considered as relevant to (e.g., as a basis for or as part of) Computer Assisted 
Diagnosis (CAD), each aspect of the invention claimed herein includes such relevancy to 
CAD apparatus, systems and methods. 

10 Section 2: Nodule Size and Shape — As a means for explanation and analysis, a 

basic model for the appearance of a small solid pulmonary nodule is defined as roughly 
spherical, with a diameter less than 1 cm. The nodule has a density (CT attenuation) 
significantly higher than that of the surroimding limg parenchyma. The nodule form may 
be confounded by other neighboring/attached structures, including vessels and the pleural 

1 5 surface. As the concem is the analytic characterization of the nodule itself, an important 
task is image segmentation the delineation of the nodule boundary and removal of other 
structures firom consideration. The image segmentation process is one of the most 
challenging aspects of any computer vision task. Segmentation methods for pulmonary 
nodules in CT scans are discussed herein. Segmentation can also be viewed as 

20 preprocessing of the image region of interest (ROI) to yield a three-dimensional 

representation of the nodule which conforms to our basic model of a pulmonary nodule. 

Consider this basic model in the case where the nodule is well circumscribed. An 
image of a well-circimiscribed pulmonary nodule is shown in Figure 9. In this case, the 
segmentation may be accomplished usiag a simple gray-level thresholding operation. 

25 Selection of the appropriate threshold involves the following issues. First, a threshold 
can be selected that best separates the nodule (foregroxmd) firom the lung parenchyma 
(backgroimd), which in the well-circumscribed case would have a bimodal intensity 
histogram. This separation may be achieved by selection of a threshold that falls at the 
midpoint between the peaks of each mode of the histogram, or by using a process (and 
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related algorithmic relationship) that considers the area under the histogram in each 
region. Thresholding methods are discussed in Section 4. 

The second consideration when selecting a threshold is in evaluating its change in 
size over time. Thus a single threshold should be selected which is appropriate across 
5 multiple studies of the same nodule. In addition, a standard threshold may be determined 
for a general class of well-circumscribed nodules such that the size of nodules can be 
compared across time and multiple subjects. 

Confounding Structures - Well-circumscribed pulmonary nodiiles with their 
characteristically high contrast to the surrounding limg parenchyma are unfortunately not 
10 the only class of pulmonary nodules commonly seen in CT studies. Measurement of the 
size and shape of pulmonary nodules is frequently confounded by the presence of other 
structures nearby or adjoining the nodule. A simple two-level model based on CT 
attenuation is insufficient to completely separate the pulmonary nodule under r- 
consideration j&om these other structures, as they frequently exhibit similar density 
15 characteristics. 

Vasculature - Pulmonary nodules, as other cellular structures in the body, require 
oxygen and nutrients to sustain growth. These are supplied to the nodule by the local 
vasculature (blood vessels). The radiographic appearance of these vessels ranges from 
undetectable in peripheral lesions, to significant in more central ones. Figure 10 shows 
20 an image of a small pulmonary nodule with surrounding and attached vessels. 

A variety of physiological factors underly the presence of vascular attachments to 
pulmonary nodules. In the context of this work, however, it is sufficient to state that 
vessels may appear nearby or significantly attached to nodules, which complicates the 
nodule segmentation problem. Method (and related algorithms) for the subtraction of 
25 blood vessels during im^e segmentation are discussed in Section 4. 

Pleural Surface - Peripheral pulmonary nodules often exhibit some degree of 
attachment to the pleural surfece (on the periphery of the lung, compressed against the 
external boundary of the thorax). One possible type of attachment is the so-called 
"pleural tail sign." This appears as a small, thin structure of nodule-like attenuation 
30 coxmecting the nodule to the pleural surface. A nodule with a pleural tail is shown in 
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Figure 1 1 . In such cases, segmentation is relatively straightforward and methods similar 
to those used to remove confounding vessels may be employed. 

A second class of pleurally-attached nodules are those that are juxtapleural. 
These nodules share a significant amount of their surface with the pleura. For this reason, 
5 delineation of the boundary between pleura and nodule can be quite challenging. A small 
juxtapleural nodule is shown in Figure 12. Segmentation techniques for juxtapleural 
nodules are described in Section 4. 

Nodule Models - Segmentation and measurement of pulmonary nodules in CT 
data benefits from a more formal mathematical model and description of each nodule 
10 type we wish to describe. The four basic solid nodule types are as follows: well- 
circumscribed, vascularized, pleural-tail, and juxtapleural. Having given a qualitative 
description of each one of these types, we may now formalize mathematical models for 
each. 

Well-Circumscribed - A well-circumscribed nodule is one that is located 
1 5 centrally in the lung, without significant connections to vasculature. A simple model for 
the nodule woxild be based on the fact that there is contrast in x-ray attenuation and, 
therefore, voxel intensity, between the nodule and the surroimdmg lung parenchyma. 
This model a set of voxels in the nodule would then include all voxels greater than or 
equal to a given intensity threshold: 

20 N= {v(x, y, z) i v{x, y,z)>T\ (2.1) 

This model would be very sensitive to noise, however, as any voxel in the region 
of interest meeting the threshold criterion would be included. A better model is to include 
connectivity between the voxels in the nodule. For each voxel in the set, then, there 
would have to be at least one neighboring voxel also in the set. Note that we intentionally 
25 exclude the trivial case where the nodule would be a single voxel. This connected nodule 
model coiild be defined in terms of the the set of voxels adjacent to a given voxel. This 
set will be denoted as adj(v(x, y, zj). Voxel coimectivity will be described below in detail. 
Our model for the well-circumscribed nodule is then the set of connected voxels greater 
than or equal to a given intensity threshold: 

39 



wo 01/78005 



PCTAJSOl/11820 



K = \y{x,y,z)\iv{x,y,z)^T)Ai3n\n^ adj(y{x,y,z)\n e Nj) (2.2) 

There is the possibility that there may be more than one coimected component A 
in the region of interest (ROI). Considering the fact that there may be more than one set 
of coimected voxels meeting the threshold criterion, we will elect to choose the largest of 
5 these sets, Nmax^ (that which has the largest volume): 

Nmaxy^rnsoiAi (2.3) 

where each of the eligible sets is dejQned as 

A = {v(x, z) I (v(x, J/, z) > r) A (3 « I « e adj(y{x, y, z)\ n^A,)] (2,4) 

Thus, we will assume that the nodule will be the largest connected component 
10 (meeting the threshold criterion) in the ROI. 

Vascularized - A vascularized nodule is one that is located centrally in the lung, 
but has significant V£iscularization (connections to neighboring vessels). Clearly, the use 
of the well-circumscribed nodule model for the segmentation and analysis of 
vascularized nodules is inappropriate, as that model of connectivity would include the 
15 vascular components as part of the nodule. Instead, we define a model which 
distinguishes nodule volume from connected vessels geometrically. 

Consider a model in wliich the vascularized nodule is the union of spheres of 
diameter 

N,^[jsix,y,z)\d,^p. (2.5) 

20 Each of the voxels in each one of these spheres must satisfy a threshold criterion, 

as W£is true for the well-circumscribed nodule model, as well as have no more than /x/2 
distance from the sphere's center. The following equation describes such a sphere, s(x, y, 
z)y centered at voxel (x, y, z)y where D is the Euclidean distance function between two 
voxel locations. 
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We will additionally define a maximuna diameter, K of any vessels which may be 
connected to the nodule volimie. Then, if we choose fi such that fi>Kwe ensure that the 
majority of each vessel volume will not be included in the nodule volume. 

5 This model of a vascularized nodule is illustrated in two dimensions in Figure 13. 

The illustration on the left shows all of the nodule and vessel area (volume in 3D) which 
meets the thresholding criterion. The illustration on the right shows how the nodule 
volume can be described as the union of circles (spheres in 3D) of diameter ja, while 
effectively ignoring the majority of the vascular component The model of a vascularized 

10 nodule as a union of translated spheres is closely allied with the notion of the opening 
operation in mathematical morphology. This topic and it's practical relation to tiie 
vascularized nodule model will be discussed below in greater detail. One of the issues of 
interest is the overestimation of nodule volume at the pohit of vessel attachment, as 
illustrated in Figure 14. Note how the spherical "structuring element" allows the 

15 inclusion of a small amount of the vessel volimie (outside the spherical nodule volume). 
For an idealized model of a spherical nodule with a vascular attachment of diameter X, it 
is possible to formally state the extent of this volume overestimation. The volume in the 
extraneous region, Vext, is the difference between the portion (spherical cap), Vfi, of the 

A 

Structuring element that extends beyond the line of length X, and the spherical cap, V^, of 
20 the spherical nodule volvime beyond the same line. These two regions are illustrated (in 
2D) in the right-hand part of Figure 14. The extraneous region, Vext is tiiat between the 
dotted line (the bovindary of the large sphere) and the external .boundary. The largest 
perpendicular distance in this cup-shaped volume is shown as v. The volume Fju is the 
region between the straight line of length A and the external boundary, while the volume 
25 V is that between the straight line and the dotted boundary of the large sphere. 

s 

Therefore, the extraneous, cup-shaped, volume is simply their difference: 

Vc-V^-t (2.7) 



41 



wo 01/78005 



PCTAJSOl/11820 



Figure 1 5 illustrates the volxime of a spherical cap. The volume as a function of r, 
. the radius of the sphere, and a, the radius of the intersecting disk is 

K(r,a)=f(r-Vr'-a')' (2r- 4^^^^) (2.8) 

Using Equation 2.8, we may substitute the appropriate values of r and a to derive 
5 the volumes of F/i , and , and subsequentiy, Vext- For the structuring sphere, of radius 
)U/2, the volume of the spherical cap is 

^^[lxl2-4it^l2f-XI2f) (p-V«/2)'-(^/2)^) (2-9) 
and for the large sphere of radius r^, tiie volume of the corresponding spherical cap is 

Vs =f^s -ylr^-^/2rJ (2r, -^r^-iXJlf) (2.10) 

The volume of the extraneous cup-shaped region, therefore, can be determined using 

Equations 2.7, 2.9, and 2.10. Note that these expressions are general for any sized 

spherical nodules, structuring spheres, and vessel diameters. More specifically, for a 

particular spherical nodule size and structuring sphere, tiie contribution of multiple 

vessels of different diameters can be determined. An example is illustrated in Figure 16, 

where a spherical nodule of 8 mm in diameter (r = 4) has a vascular attachment of 

varying diameter (X). Each of the curves in the graph corresponds to a A, value of 0.75, 

I.O, or L25 mm. The graph illustrates the percent overestimation in nodule voliune when 

the diameter of the structuring kernel used (jiJ) is varied from 1.25001 to 2.5 mm. In the 

worst case (X = 1.25 mm, fx = 1.25001 mm), where the structuring kernel is barely large 

enough to detach the vessel, the volume overestimation is less than 0.19%. As the kemel 

diameter is increased, this value decreases logarithmically. 

Pleural Tail - The third variety of nodule is near the pleural surface and 

coimected to it by a thin structure, called a "pleural tail." While it may be advantageous 

in some respects to treat this group differently than those in all other groups, it may be 

modeled as one of the other two groups, as a vascularized nodule, or as a juxtapleural 
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nodxile. In the former case, consider a model similar to that of the vascularized nodule^ 
where the nodule region is the union of spheres of diameter wherein each voxel meets 
a fundamental threshold criterion 

N,^[jsix,y,z)\d,^tx (2.11) 

5 This expression is sufficient to exclude all of the pleural tail except for a small, cup- 
shaped region. A fiill discussion was given in the context of vascularized nodules set 
forth above. When the volume of the pleural tail is to be considered jointly with the rest 
of the nodule, the juxtapleural nodule model tray be used, described in the next section. 

Juxtapleurai - In juxtapleural nodules, a significant amount of the nodule 

10 periphery is connected to the pleural surface. Figure 17 illustrates the basic model. Four 
regions are defined: N, the nodule region; LP, the surrounding lung parenchyma; PS, the 
pleural surface; and TW, the thoracic wall region. In CT images, the nodule region and 
the thoracic wall region both exhibit high-attenuation characteristics. Therefore, while 
simple thresholding is able to separate the nodule j&om the lung parenchyma, it is 

15 insufficient to separate the nodule from the outer thoracic region. In Figure 17, a dashed 
line between points a and b has been included to indicate one possible boundary between 
nodule and thoracic wall regions. Physiologically, the pleural surface separates the lung 
volume (containing the parenchyma) and the thoracic wall. When a juxtapleural nodule 
is present, however, there may or may not be invasion of the pleura. Whether the nodule 

20 has invaded the pleura or is merely adjacent to it, voxel density (intensity) information is 
rarely sufficient to distinguish these regions. 

In our model, we would like to consider the nodule region to be only that volimae 
which is unlikely to be part of the thoracic wall. This can be defined geometrically. We 
may define points a and b in the two-dimensional cross-section shown in Figure 17 to be 

25 those points of maximum concavity for a given scale. The notion of scale in the 
estimation of concavity here is quite important, as the pleural surface boundary is likely 
to have local concavities due to image sampling geometry and also physiologic/anatomic 
variation. Techniques for determination of the boundary between the nodule region and 
thoracic wall will be discussed below. 
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Section 3: Data Acquisition «- Data acquisition and management are non-trivial 
aspects to lung cancer research using computed tomography. All CT data used in this 
work were collected as part of the IRJB-approved Cornell Early Lung Cancer Action 
Project (ELCAP), The data were acquired using GE Highspeed and LightSpeed CT 
5 scanners using the protocols described the Background and shown in Figure 3 . 

Each two-dimensional CT image is a 512 x 512 hnage with 12-bit pixels, 
providing a dynamic range of 4096 gray-levels. The pixels are backed into two-byte 
words (16 bits), simplifying the data alignment for storage and processing, A given CT 
examination may generate as many as 200 of these images, given several series and 
10 associated reconstructions. A quick calculation illustrates the enormity (by current 
technological standards) of the storage problem. 

Sizeintage = 512 •512-2 bytes = 512 KB (3.1) 
Sizeexcw, = 200 • 512 KB = 100 MB (3.2) 
While the raw data for most studies is on the order of 50 MB, additional storage 
1 5 must be available for processed images, results of image analyses, and special 

visualizations. Given the large data storage and management reqxnrements of this 
project, standardized mechanisms were developed for image accrual, processing, and 
archiving. 

DICOM - The Digital Imaging and Communications in Medicine (DICOM) 3.0 
20 [13] communication protocol and image storage standards can be used as the interface 
to the system. These are well-defined standards which allow the acquisition of data 
from a variety of sources, providiag a system and vendor independent interface to the 
analysis system. 

VisionX - The VisionX system, developed at Cornell University, is a collection 
25 of computer vision and image processing libraries which have been used in a wide 

range of image analysis research. The VisionX development system was the basis for 
the method, system, and software developed in this invention. 

Software for the transmission of DICOM images to and from PACS (Picture 
Archiving and Communication Systems) systems and conversion to and VisionX format 
30 was developed. In VisionX format, the DICOM header information, selectively blinded 
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with respect to patient-identifyiiig information, is retained to facilitate retransmission of 
data and/or analyses to standard PACS systems. 

Isotropic Resampling -Most CT data are sampled in an anisotropic space, 
where the resolution in each of the three dimensions is not equal. The in-plane (^—y) 
5 resolution on modem scanners is typically better than 1 .0 mm in both dimensions. The 
axial (z) resolution, however, may be anywhere from near 1 nam to as high as 20 nam. 
A typical high resolution scan protocol provides 0.5 mm in-plane resolution with 1 mm 
axial reconstructions, produciug voxels that are (0.5 x 0.5 x 1 .0 nun). This poses 
several problems. First, although high-resolution 2D measurements can be made in any 

1 0 one of the highly-resolved CT slices, 3D measurements are hindered by the less- 
resolved axial dimension. Special anisotropic image filters need to take into account the 
different axial resolution, which becomes computationally cumbersome. A second, 
more important problem to consider is that of partial volume effects. 

Partial Volume Efi^ects - The partial volxmie problem describes the quantized, 

15 discrete spatial nature of the CT image. Consider the following two-dimensional 

example. A unit circle (radius = 1) is perfectly sampled on a regular grid. Pixels are set 
if a sufficient percentage of their area corresponds to the interior of the circle. We may 
assume this threshold to be 50% of a pixel's area. The resolution of the system, at the 
coarsest setting is equal to the radius of the circle. This is illustrated graphically in 

20 Figiire 1 8. For simplicity, most of the following discussion will concern only one 

segment of the circle, the upper left 90"^ segment, and the corresponding quadrant of the 
square in which the circle is inscribed. The quadrant may be divided at any number of 
regular intervals to dejQne an appropriate sampling grid, but for this example, we will 
choose this division factor to be powers of two for the number of divisions in the x and 

25 y dimensions, with n the appropriate exponent of 2. In Figure 1 8 , for example n is 
equal to zero, and the quadrant is divided into one large pixel, with area equal to 1 . 
Figure 19 illustrates the cases where n is equal to 1 and 2. 

With this model in place, we may discuss the accuracy of area measurement of 
the circular segment as a function of sampling interval (image resolution). The true In 

30 the area of tiie circular segment is 
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— 0.7854 
4 

In the case where n—l, the area is estimated to be 1 .0, as more than 50% of tiie pixel is 
covered by the circle. When n is increased to 1 and 2, the area estimates improve. In 
the first case, three of four pixels are set, yielding an area estimate of 0.75. For « = 2, 
5 1 3 of the 1 6 pixels are sufficiently filled by the circle, giving an estimate of 0.8 1 25. 

Figure 20 illustrates the cases in which n is equal to 1 , 2, 3, and 4, showing only 
the quadrant of interest. In each case, the pixels that would be considered members of 
the circle are shaded. It is obvious firom the figures that as the degree of subdivision of 
the sampling grid increases, the error in area estimation of the circle decreases. It is 
10 also possible to produce boimds for this error as a function of the degree of subdivision. 
Consider, however, the case where the unified quadrant is the original pixel size in the 
acquired image. We may then cast our discussion in terms of the degree of resamplings . 
of the image space. We may define the supersampling ratio, to be the number of 
divisions of the original pixel size in each of the x and y dimensions. The n-th term of 

15 the series describing then, is simply 

% 

Sn^Z"" (3.3) 

Similarly, the number of pixels,/?, in the quadrant is given by 

Pn^2'^ (3.4) 

20 Now, consider that the error in area estimation is the sum of errors made in 

representing the circle boundary, and that we therefore only need to consider the pixels 

through which this boundary passes to evaluate the error. Each of these pixels will 

each contribute to either an overestimate or underestimate of the true area. For 

example, in the upper left illustration (n-l) \sx Figure 20 , the upper right and lower 

25 left pixels are set, each overestimating the true area of the segment. The upper left 

pixel however, is unset, contributing to a local underestimation of the area. Again, it 

should be clear that only those pixels through which the true boundary passes 

contribute to this error. Pixels that are definitively in the image foregroxmd (set) or 

backgroxmd (unset) will be correct. In the upper right illustration (n = 2) in Figure 20 , 
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there are 8 definitive foreground pixels and one definitive background pixel. The 
remaining seven, those on the boundary, each contribute some positive or negative 
error to the overall area estimate. 

It is possible to formally define the number of such boundary pixels in this 
model, those for which the relative proportion of foreground will determine whether 
they are set. The number of these boxmdary pixels, b, in this quadrant is defined as 



The first few terms in the series, { 1, 3, 7, 1 5, 3 1 }, corresponding to 0 < « < 4, describe 
the numbers of upper left quadrant boundary pixels shown in Figures 18 and 20. For 
example, in the lower left illustration in Figure 20 (n = 3), the supersampling ratio, is 
8, the nimiber of pixels, j7, is 64, and the number of boundary pixels^ b, is 15. 

The magnitude of the contribution each boxmdary pixel to the overall error is 
boimded by 0.5 times the area of the pixel, as we chose 50% coverage as the threshold 
for determining whether a pixel is set. Therefore, a conservative upper bound on the 

error for estimatmg the area of the circxilar segment, € can be expressed as 



where a is the area of a pixel. Given that this is a unit circle, and therefore the quadrant 
has unit area, the area of each pixel is 




(3.5) 



2 



(3.6) 




(3.7) 



Thus, the error bound for a given value of n may be restated as the series where 




(3.8) 



and, after some manipulation. 




(3.9) 
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Each of the two tenns in this series converges to zero, with the latter converging much 
faster* Essentially, as n increases by one, the number of boundary pixels basically 
doubles, while their size, and thus their error contribution, is divided by four. This 
leads to an error bound that is approximately halved with each increase in n. 
5 With this two-dimensional example, we can see how perfect interpolation, or resampling, 
of the image space can reduce errors in size estimation. In three dimensions, isotropic 
resampling (supersampling to an isotropic space) allows segmentation decisions to be 
made on a super-resolved, or supersampled, grid, allowing more accurate, consistent 
boxmdary decisions to be made. The intensity in each image voxel is interpolated to 

10 estimate the intensities in each of the sub-voxels in the supersampled space. Thus, each 
original voxel iatensity value is responsible for several in the new image space, 
mitigating partial volimie effects ia the original data, as a more precise spatial location of 
the desired gray-level transition (boundary) can be identified. 

Supersampling - Three-dimensional supersampling is the process by which the - 

15 image is given a new representation in units bearing a higher spatial frequency than 
those in the original image, resulting ia an improvement in image resolution. The 
original image may be anisotropic or isotropic and supersampled to either a more 
refined anisotropic or isotropic space. The process involves computing new voxel 
values at every location in the new, more refined coordinate space. Each of these 

20 values may be computed using one of a variety of interpolation schemes including 
trilinear, tricubic, and B-spline intetpolation. 

Trilinear interpolation - In the present invention, supersampling was conducted 

♦ 

using trilinear interpolation of each voxel value at location yn, Zn), where each new 
voxel value can be defined in terms of the voxels' eight nearest neighbors in the original 
25 space, v(jc„ yj, z0, where 

i=[0|l] 

7 = [0|1] (3.10) 
^=[0|1] 

For notational simplicity, we define three distances between v(xn, yn,, Zf) and v(x„ yp zp) 
30 in terms of the perpendicular distance along each of the jc, j/, and z Cartesian axes. 
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dXi= \xn-Xi\ 

(fyj=' \yn-yj\ (3.11) 

We also define the resolution in each dimension in teims of the distance between 
5 adjacent voxels in each dimension. 

ry= \yi-yo\ (3-i2) 

rz= Jz;-xo| 

Given the above, the value of the new voxel based on trilinear interpolation 
10 is 

+ v(x„>;o>^o) • dx^lrx-il-dy^/ry) • (1 - dz^/rz) 
+ y(Xa,yi»Za') • (1 - dxjrx-dy^/ry ■ (1 - dz^lrz) 
+ v(Xo, 3/0,21) 'Q-dxo/rx)- il-dyo/ry) • dz^/rz) 
15 +vixi,yo,z^)'cixjrx-0.-dyo/ry)'dzo/rz) 

+ v(Xo ,yi,zi) (l-dxjrx'(l-dy^ fry) Q-dz^ jrz) 
+ viXi,y^,Zo)'dxJrxdyo/Ty'il-dzQ/rz) 
+ v(xj , , Zi) • i&o /rx -dy^lry'dz^l rz) (3.13) 
An alternative, more concise expression for the value of v(x«, y^ Zw) would be 

20 vix„,y„,z„) = Y,i^T,C^-dxJrx)il-dyj/ry)(l-dz,/rz) vix,,yj,z,) (3.14) 

/aO y»o JtM) 

Tricubic interpolation - Supersampling was also conducted using tricubic 
interpolation which generates a new voxel value at location y„y Zk,) based on the 64 
nearest neighbors in the original space, v(Xi, yi, Zk\ where Uj\ and k can have the values - 
1, 0,1, 2. Briefly, the new value can be expressed as 

25 v(x„.;.„,zJ = ttti'(^^)^'(V^>^<T^)^<*''^^'^*> (^-l^) 
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where P(f) is a cubic polynomial interpolation function with the following 

specification: 

'l-2(/2)+|/'| : 0<|/|<1.0 
4-8-|r| + 5.r'-|r^| : 1<|^|<2.0 (3.16) 
0: 1^1 > 2.0 



Although, the inventors have used the methods set forth hereinabove, the present 
5 invention includes any method of supersampling which provides sufficiently 
improved image resolution. 

Section 4: Segmentation - Image segmentation is the process by which an image 
is divided into regions corresponding to objects and backgroimd. In the analysis of 
pulmonary nodules in radiographic images, segmentation involves the differentiation of 
10 nodule from non-nodule regions. This chapter details algorithmic solutions to the nodule 
segmentation problem, both for well-circumscribed nodules and for those with 
attachments to vasculature or the pleural surface. 

Basic Segmentation - The automated (and manual) analysis of pulmonary 
nodules is intimately tied to the image segmentation problem. In order to take 
15 measurements of nodule size, shape, and density characteristics, one must first determine 
which image regions belong to the nodule and which do not. A first order solution to this 
problem can be described using a two-level model. 

Two-Level Model - Considering the region of interest (ROI) as composed of two 
classes of voxels, those that are part of the nodule (foreground) and those that are not 
20 (backgrotmd). Figure 21 shows a single 2D image of a nodule. The bright, high- 
attenuation region is the nodule, surrounded by a darker, low-attenuation region 
corresponding to lung parenchymal tissue. 

Examining the frequency distribution of gray-level intensities in this image, it is 
found that they exhibit a bimodal distribution. A histogram of the image gray-levels (in 
25 HU) is shown in Figure 22. The large mode centered at around -850 HU represents the 
low-intensity pixels ia the background. The smaller mode centered near 0 HU 
corresponds to the higher-intensity nodule region. 
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Consider the two-dimensional image function,/ (x, y). A simple two-level model 
can be defined where foreground and background pixel classes may be separated by 
intensity. A transformation r maps pixels in two one of two classes, 0 or 1 . A simple 
transformation that relies on the two classes having complete separation by intensity can 
be described as follows: 



This gray-level intensity transformation is commonly known as thresholding. The 
challenge in separating the classes in this two-level model, then, is the choice of an 
appropriate threshold value, Several methods of threshold selection are described in 
1 0 the next section. 

Thresholding - The goal of threshold selection is to determine a threshold which 
optimally separates backgroxmd from foregroxmd regions. While this may seem 
straightforward for a simple image (such as that shown in Figure 21), the optimality 
criteria may change in a variety of circimistances. Here are some examples: (a) 
1 5 maximize separation of two modes in a bimodal histogram derived from a single image; 

(b) maximize separation of two modes in a bimodal histogram of an entire 3D ROI; 

(c) maintain a consistent segmentation method between scans of the same nodule over 
time; and (d) maintain a consistent segmentation method across nodules in multiple 
subjects. 

w 

20 The threshold selection problem for a single 2D nodule image may be solved 

using histogram analysis and selection of a threshold vedue which best separates the two 
modes of the (e3q)ected) bimodal distribution. Many mathematical relationships (e.g., 
algorithms) have been proposed for the global thresholding problem (a single threshold 
for the entire image), as well as for local thresholding (thresholds depend on local image 

25 properties in addition to pixel intensity) [24, 36]. A simple global threshold selection 
algorithm will now be described, based on the method of Ridler and Calvard [68]. 

Algorithm 4,1 (Iterative Threshold Selection) 

Select a candidate threshold value T for histogram H(p) 
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while ((Am, >e)i|(AiU2 >€)) 

i?i =Hip):0<p<T 
^Hip):T<p<co 

Mip = Ml 
M,=(Si?0/KI 

iLi2=(2i?2)/|^2| 

r-(Mi+M2)/2 
AMi=iUi-i"ip 

AM2 =A^2-M2p 

end 

Algorithm 4,1 iteratively selects a threshold based on the mean intensities (juu 
JU2) of the two histogram regions (Ru ^2) divided by T. At each iteration the new value 
of T is computed as the mean of the mean values of each region. 

r = (jUi+M2)/2 (4.2) 

The termination condition is that the means of each region change by no more 

than €, which is somewhat application dependent. The result of Algorithm 4.1 on the 

2D image shown in Figure 21 is a threshold of -476 HU. The resulting segmented 
image is shown in Figure 23. 

This method works quite well for bimodal distributions where the modes 
encompass similar area and have minimal overlap. One factor affecting the relative 
area in each of the modes, of course, is ROI selection. If the nodule is surrounded by a 
large amount of limg parenchyma (i.e. the ROI is large compared to the nodule itself), 
the high-attenuation mode may not have sxifficient area to be distinguished above the 
contributions of image noise and other high-attenuation, non-nodule structures. This 
becomes particularly important when analyzing three-dimensional regions. 
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One of the issues in threshold selection is that lower threshold values produce 
segmented nodiiles which appear larger, as the interface (margin) between nodule and 
lung parenchyma varies widely in attenuation. Figure 24 illustrates the effect of 
changing the threshold value T in 100 HU steps from -800 HU (top left) to 100 HU 
5 (bottom right) when segmenting the nodule shown in Figure 2L 

In order to make accurate estimates of nodule growth, the segmentation should 
be consistent across multiple scans taken at different times. In growth estimation, the 
absolute error in size measurement is less significant than the relative error, considering 
multiple studies of the same nodule. This is due to the fact that growth estimates are 

10 computed based on size ratios, not on absolute size. For this reason, a single threshold 
value should be used in segmentation of each of the scans, finportant caveats, however, 
are that the scanner must be well-calibrated and the scan protocol (dose and resolution 
parameters) fixed for the observed attenuation values in each scan to be comparable. 
Growth estimation is covered in detail in Section 5. 

15 It may also be beneficial to attempt a standard thresholding policy across 

nodules in different subjects. This is a more difficult problem, however, as mean 
nodule attenuation varies somewhat with nodule type. Non-solid nodiiles, known as 
ground glass opacities (GGOs), and lesions with mixed solid and GOO components 
frequently require a much lower threshold for accurate segmentation. 

20 A more general problem with global thresholding methods is that regions of 

similar intensity may be spatially disparate, and not represent components of the same 
object. For example, nearby vessels or other structures of nodule-like attenuation may 
be included in a thresholding-based segmentation. Image noise may also be of 
sufficient intensity to contribute to the output of a thresholding operation. These 

25 structures should not be included in the nodixle region to be considered for fiirfher size 
and shape analysis. 

Connected Component Analysis - One solution to the problem of disparate 

high-intensity regions contributing to the segmented nodule volume is the use of 

connected component labeling and analysis. In connected component labeling, each 

30 pixel/voxel is assigned an integer value indicating the unique 2D/3D connected 
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component of which it is a member [24]. These comiected components include all 
voxels through which a connected path can be followed. Such a path connects adjacent 
non-zero voxels that share at least one vertex. In other words, a voxel vQc, z) belongs 
to the same connected component as any of its 26 neighbors in the set iV, where 



C t 



X ax—l^x <jc+l; 
y ^y-~l<y <>'+l; 
z 3z—\<tZ 
(x\y\z)^{x,y,z) 



S and the value of the neighboring voxel is non-zero. 

The use of 3D connected component analysis allows noise and extraneous 
structures to be removed from the nodule image data. The result of selecting a single 
connected component is an object that is contiguous in 3D space. 

The first step in connected component analysis is the comected component 
10 labeling itself There two general algorithms that can be used to accomplish connected 
component labeling in three dimensions, one recursive and one iterative. 

Recursive Connected Component Labeling - The recursive connected 
component labeling algorithm takes as input a three-dimensional unage t> and produces 
a labeled image /, where each non-zero voxel value is the index of the component of 
15 which it is a member. The algorithm has two sections. The outer section loops over all 
voxels t)(x, y^ z) in the input image. When an object voxel is encoimtered (voxels with 
zero value are considered as background), the component label L is incremented and the 
recursive labeling function labelQ is called. This recursive function examines each of 
its 26 neighbors for connectivity. If a connected voxel is found, it is given the same 
20 label in a recursive call to labelQ. 

A recursive algorithm is quite elegant in that it labels each component 
completely and in tum. The primary disadvantage of such an algorithm, however is that 
each recursion places state on the system stack which is not removed until the 
termination cases are reached and the recursive calls removed firom the stack. 
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Iterative Connected Component Labeling - To eliminate the memory 
demands associated with connective component labeling, it is possible to formulate an 
iterative algorithm, albeit at the cost of complexity and execution time. 

The iterative method traverses the image in order, in increasing x, y, and z. If 
5 connectivity is estabUshed between the current voxel v(x, z) and any of the 1 3 

previously searched (o (X'l..,x + l,y'l...y+l,z-V) (9), o (x - 1 . . .x + l,j;-l, 
z) (3), u (x - 1, y^ z) (1)), all possible connections are recorded in a table. At the end of 
the iterative search through the image space, this table must be reconciled so that 
coimected components originally assigned different labels that are subsequently found 
10 to be connected are given a consistent label. 

Size and Spatial Analysis - Once connected component labeling has been 
performed, more complete connected component analysis can take place. In 
segmentation of pulmonary nodules, several selection criteria are commonly used in 
connected component analysis. These criteria are used to (i) select the component of 
15 greatest voliune, (ii) select all components greater than or equal to a specified volume 
threshold, and (iii) discard components within a specified distance of the ROI 
boundaries. 

These three criteria are used in the following way. 1) The object of greatest 
volume in the ROI is typically the nodxile. 2) In some cases, more than one object of 
20 high relative volume in the ROI may need to be selected. 3) Lastly, segmentation of 
nodules near other large structures (e.g. pleural surface) may not be the object of 
greatest volume. In these cases, the extraneous object is typically close to (or adjoining) 
the ROI boxmdary. 

Thus, image thresholding and connected component analysis can be used for 
25 segmentation of pulmonary nodules when they are of relatively high intensity compared 
with the surrounding lung parenchyma. Such methods are insufficient when nodules 
are attached to other structures of similar intensily (e.g. blood vessels). 

More advanced segmentation techniques are required for isolating the nodule 
volume in these cases. Such techniques have been developed in this invention and are 
30 based on morphological filtering. 
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Morphological Filtering -Mathematical morphology is the study of non-linear 
image filters based on geometric constraints. This field was first developed by G. 
Matheron and J. Serra in the late 1960's and popularized in computer vision and image 
processing in the subsequent decades [74, 75]. Morphological operations are useful in 
edge-detection, shape-detection, skeletonization, and other image segmentation tasks. 
The fimdamental morphological operations, erosion and dilation are based on 
comparison between small regions of an image and a structuring element or kernel, 
translated across the image in much the same way as the convolution kernel used in 
many other filtering operations. Morphological operations were first defined to operate 
on binary images, and whereas they have been extended to gray-level morphology, the 
c\zrrent discussion concerns the binary morphological transformations [74, 35]. 

Erosion is the set of all voxels in an image A that can include a small kernel, B. 
It can be expressed as 

AGB^ {v(x,y.z) B c= a} (4,3) 

where B(x^^) indicates the kernel B centered at voxel (x, z). 

Dilation is the set of points where the translated kernel B(x, y ,z) would touch 
the image ^. This can also be expressed as 

A®B^ {v(x, z) 3 } n A ^ 0 (4.4) 

An additional operation of interest is reflection, 

A' = {v(-x-y-z)\v(x, y, z) e a) (4.5) 

which is geometric reflection of all voxels about the center of a structuring element or a 
complete image. 

Two additional morphological operations are defined as each composition of the 
previous two, as erosion and dilation do not form an inverse pair. Opening, is the 
erosion of an image ^4 by a kernel B followed by dilation of the result using the same 
kernel: 
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Aj, =iAeB) ®B (4.6) 

Sometimes referred to as a "rolling-ball" filter, the result of a morphological 
opening operation can be viewed as the interior region of an object covered by 
translation of the kernel into all areas where it remains completely enclosed. An 
5 opening operation has the effect of breaking small connections between regions and 

removing sharp edges. 

The dual operation^ closing, is a dilation followed by an erosion: 

^iA®B)eB (4.7) 

* 

Closing can be viewed as an opening operation performed on the background 

1 0 (non-object region) of an image using the geometric reflection of the structuring kernel, 
B\ Graphical examples of morphological erosion, dilation, and opening are shown in 
Figure 25. Vote the effect of the structuring element on each of the surface features in 
the dilation example. An opening operation gives the resultant surface the 
characteristics of the structuring kernel; in this case, the circular characteristics of the 

1 5 disk-shaped kernel. In addition, note that to remove thin structures cormected to a 

surface (e.g. vessels), the kernel must be. of a size greater than that of the structure to be 
removed. As suggested by this example, in the segmentation of puhnonary nodules, a 
morphological opening operation may be used to reduce the likelihood that noise 
artifacts and small vessels be considered part of a nodule. 

20 Vascular Subtraction - Global methods for removing the puhnonary 

vasculature, based on three-dimensional region growing, tree-traversal algorithms, and 
other techniques have been described [90, 16, 79]. The goal of these global methodsis 
to simplify the nodule detection problem. For the characterization of pulmonary 
nodules using high-resolution data, however, local methods of vascular subtraction are 

25 more appropriate. Since only a small subset of the lung volume is normally acquired in 
a high-resolution diagnostic study, global methods tracing the vascular branches from 
the pulmonary artery and veins onward are impossible. Even if high-resolution data of 
the entire lung volume were available, this approach would likely be impractical. 
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Furthermore, global or semi-global region growing schemes based on voxel intensity 
alone have the danger of removing small nodules when they exhibit vascular 
connections. 

Methods - As an alternative to global region growing techniques, we use a local 
filtering approach based on mathematical morphology. This method is based on the 
model described in Section 2. The initial morphological processing in our segmentation 
algorithm consists of an opening operation based on a spherical kernel. This kernel is 
passed over the input data in a convolution-like filter for both the erosion and dilation 
steps, followed by connected component analysis to retain the nodule volvmie and 
discard vascular components that were initially disjoint or were disconnected via the 
opening operation. 

One disadvantage of morphological opening in this application is that it has a 
"smoothing" effect on the nodule surface, as surface details are removed, such as fine 
spiculations. Consider the illustration in Figvire 26. The illustration on the left depicts a 
pulmonary nodule with a vascular attachment. In the center illustration, the translation 
of a circular kernel is shown, depicting a 2D morphological opening operation ("rolling- 
ball filter"), hx the result shown on the right, notice that while the vessel has been 
removed (the desired effect), some of the surface features present in the original nodule 
have been deleted, as they were smaller than the structuring kemel. 

While the need to remove vessels firom consideration is important, we would 
prefer not to smooth away the very nodule svirface characteristics we hope to analyze. 
To compensate for this smoothing effect, we may perform an iterative constrained 
dilation (morphological opening followed by geodesic dilation) process to "regrow" 
these features. The entire morphological filtering process is as follows: 

Algorithm 4,2 (Iterative Morphological Filtering) 

Begin with an initial binary image / 

J - (le Sd) ® Sd {Perform opening using a spherical kemel Sd of diameter d\ 

Perform connected component analysis, keeping the component of interest 
while s > = 2 {Nxmiber of useful dilations} 
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J= J © & {Perform dilation using a spherical kernel of diameter s} 
J=JAI {Perform a voxel-by-voxel logical AND} 



end 



5 



This technique restores the detailed surface features of the nodule without 



regrowing vessels more than L — d voxels from the surface. In the first iteration, the 
dilation with a sphere of diameter d extends the surface at most d/2 voxels in the 
direction of the vessel. In each subsequent iteration, the surface may be extended by, at 



most, half the distance of the previous iteration. The upper bound on the growth of each 
10 vessel, therefore, is a distance of d voxels, as can be seen from the following geometric 



In addition, the logical AND operation guarantees that all the features that are 
grown were present in the initially thresholded image. Figure 27 illustrates the vascxjlar 

15 subtraction using morphological opening and connected component analysis, but without 
iterative filtering to regrow surface features. Two 3D shaded surface models are shown. 
The left image shows the result of a basic segmentation method based on thresholding. 
The right image shows the result of a 3D morphological opening using a spherical kernel 
of an appropriate size to remove the vessels connected to the nodule (as well as others in 

20 the ROI), via conoaected component analysis. Note that the surface of the segmented 
nodule is significantly smoothed when compared with the original thresholded data. 

In comparison. Figure 28 illustrates the result of vascular subtraction via 
Algorithm 4.2, which adds the iterative morphological operations needed to regrow the 
nodule surface features. Note that the surface features present in the original thresholded 

25 data have been restored. 

Figure 29 shows another example of vascular subtraction via iterative 
morphological filtering. In this example, the major vascular component is removed, in 
addition to several minor vessels (note the small vessels protruding from the top surface 



* 



senes. 




(4.8) 
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of the threshold-segmented nodule (left). Again, the vascular subtraction was achieved 
without significant effect on desirable nodule sur&ce features. 

Sensitivity Analysis - The two critical choices in the segmentation algorithm for 
vascularized nodules are the gray-level threshold and the structuring kemel diameter, 
5 Methods for threshold selection were discussed hereinabove. In this section, we will 
examine the sensitivity of the morphological filtering methods used to remove vascular 
structures from the thresholded nodule ROIs. 

The diameter of the structuring kemel used in Algorithm 4.2 affects both the 
initial morphological opening of the scene, but also the sizes of kernels used to regrow 

10 the surface features. The most important consideration when choosing the kemel size is 
that the sizes of nodules and vessels vary, as described in our vascularized nodule model 
(Section 2). In particular, the cross-sectional diameter of attached vessels (X firom our 
model), may vary considerably firom case to case. Still, although we may choose a kemel 
of diameter rf, such that it is likely to be larger than most vessels, overestimation of the 

1 5 appropriate kemel size may lead to overestimation of the nodule voliome at those points 
where considerably smaller vessels were attached. This relationship was expressed 
mathematically in Equation 4.8, where the upper boimd on "regrowth" of a vessel is 
bounded by the diameter of the stmcturing kemel, d, A complete expression for this 
overestimate was derived for the spherical vascularized nodule model in Equations 2.7, 

20 2.9, and 2.10. Here, we will illustrate graphically for a three-dimensional synthetic 
model. In this example, a synthetic nodule was constructed having a diameter of 25 
voxels, witii an attached vessel of 5 voxels in diameter. Given our standard 0.25 mm 
isotropically resampled data, this would correspond to a nodule diameter of 6.25 mm and 
vessel diameter of 1 .25 mm, well within the typical range in this study. 

25 Using this synthetic nodule, the diameter of the stmcturing element used in 

Algorithm 4.2 was varied and the resultant segmentations evaluated visually and 
numerically. Given the characteristics of the iterative vascular subtraction method, there 
are three regions in which values of this parameter may fall: (i) d is too small - the vessel 
is not removed; (ii) is in an acceptable range - the vessel is correctly removed; and (iii) 

30 ^ is too large - the vessel is removed, but regrown to a significant degree. 
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Three examples from this experiment are shown in Figure 30. On the left (d = 4), 
the kernel is too small to remove the vessel in the initial morphological opening. In the 
center image, the kernel is in the acceptable range for a proper segmentation (in this case, 
d = 5). On the right, the kernel is somewhat too large, resulting in a significant amount of 
5 vessel "regrowth". 

An e3q)eriment was performed to assess the behavior of this segmentation 
technique using synthetic vascularized nodules (a sphere with vessels of varying 
diameters). These models were segmented using kernels of increasing size, and the 
degree of volume overestimation recorded. Figure 3 1 shows the relative volume of the 
10 segmented nodule model as a function of vessel diameter, and kernel diameter, d. 

Each of these parameters are expressed with respect to nodule spherical diameter, dg. The 
data for each vessel diameter are shown beginning with the kemel size that is sufficiently 
large to exclude the vessel. 

These data show that, although the degree of volume overestimation increases 
1 5 with kemel size (over the minimmn diameter required to exclude the vessel), the effect is 
still within 4% for kernels as large as 60% of the diameter of the nodule and for vessels 
as large as 30% the nodule. Still, the overall goal is to identify that range of kemel sizes 
that leads to good nodxjle segmentation for a wide variety of in vivo nodules. 

To this end, additional experiments were performed to test the sensitivity of the 
20 vascular subtraction technique (Algorithm 4.2) to structuring kemel diameter in 21 in 
vivo pulmonary nodules. For each nodule, the diameter of the structuring kemel was 
varied from 1 to 20 voxels. In each of the resultant segmentations, the degree to which 
the vessel or vessels (most nodules had more than one vascular attachment) were 
removed was evaluated. 
25 Figure 32 illustrates the results of these experiments. The graph shows the 

relationship between structuring kemel diameter and the percentage of nodules correctly 
segmented. Recall that, for a given nodule and set of attached vessels, there are usually a 
range of values of d that lead to a correct segmentation. 

In these results we may note that more than 80% of the in vivo vascularized 
30 nodules were correctly segmented by using a fixed stmcturing kemel diameter of 6 
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voxels (1.5 mm). The graph also exhibits a bimodal appearance. This is due to the fact 
that while the majority of vascularized nodules studied were connected to vessels smaller 
than 2 mm in diameter, there was an additional group of nodules with larger vascular 
components, resulting in the secondary peak eAd=10. Therefore, although a static 
5 structuring kemel size may be chosen that yields reasonably good results over a wide 
range of nodules, a better approach is to choose a stmcturing kemel of the appropriate 
size for a particular case. 

Pleural Surface Removal - Consider the basic model of the pulmonary nodule 
described in Sections 2. It describes two classes of pulmonary nodules that exhibit 

1 0 coimectivity to the pleural surface, those with a pleural-tail and those that are 

juxtapleural. Nodules exhibiting a pleural tail, a thin structure similar in size and density 
to the rest of the nodule, can be treated in much the same way as juxtapleural nodules. 

Juxtapleural nodules, those that share a significant amoimt of their periphery with 
the pleural surface, require a different segmentation approach. An effective segmentation 

15 of a juxtapleural nodule is one that eliminates pleural and other thoracic components 

fi-om the nodule volume. One approach to the detection of juxtapleural nodules has been 
to use 2D morphological opening with a spherical kemel within the chest wall, and 
subtraction from the lung volume [4]. 

Methods - A two stage approach is used in the segmentation of juxtapleural 

20 nodules. First, the orientation of the pleural surface in the ROI containing the nodule is 
determined using three-dimensional moment analysis [65, 64]. The method of moments 
is described in Section 6. For the purposes of the current discussion, this method is used 
to determine the angles describing the orientation of the pleural surface at the point where 
the juxtapleural nodule is attached. Once this orientation has been determined, a 

25 stmcturing kemel is generated with the appropriate size and orientation such that an 

opening operation can be used to detect the majority of the pleural surface and chest wall, 
while excluding the nodule. This kemel is disk-shaped and large enough so that it will fit 
only within the chest wall region and not the nodule portion of the image. Three- 
dimensional image subtraction is then iised to remove these external stmctures firom the 

30 original image. Lastly, the remaining pleural components not detected in the opening 
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operation are removed during the application of Algorithm 4.2. The complete method is 
described as Algorithm 4.3. 

Algorithm 43 (Pleural Surface Removal) 
Begin with an initial binary image / 
5 Using moments, determine the orientation of the pleural surface 

Generate a disk-shape kemel D, oriented parallel to the pleural surface 

J = (/ e D) © £) {Perform morphological opening using D} 
K-I- J {Perform image subtraction} 

Continue with iterative morphological filtering {Algorithm 4.2} 

10 

An example of the pleural surfiace removal techixique is shown in Figure 33, 
The nodule, approximately 5.5 mm in diameter, is removed from the pleural surface 
prior to volumetric and shape characterization. The segmentation process was as 
follows. 

1 5 The method of moments was used to compute the orientation of the pleural 

surface and attached nodule. Once the orientation of the pleural surface was 
determined, a disk-shaped kemel of 34 voxels (8.5) mm was generated and used in a 
morphological opening to identify those voxels in the pleural wall not containing the 
nodule. The pleural surface component was subsequently subtracted from the image, 

20 leaving the nodule and a small amount of pleural surface not identified by the opening 
operation. The resviltant image was then segmented using iterative morphological 
filtering using a spherical kemel of 5 voxels (1 .25 mm) in diameter, to remove the 
remaining elements not belonging to the nodule volume. An additional example of this 
segmentation method is shown in Figure 34. As in the previous example, the nodule is 

25 automatically removed from the pleural surface. 

Sensitivity Analysis - It should be noted that the success of the pleural surface 
removal algorithm is somewhat dependent on the geometry of the nodule ROI. This is 
due to the fact that estimation of the orientation of the nodule-pleura interface is based 
on moment calculations involving the thresholded ROI. In experiments performed on 
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in vivo juxtapieural nodiiles. Algorithm 4.3 was able to perform a correct segmentation 
in 72% of the cases. When a user was able to manually specify the orientation of the 
disk, the results improved, allowing for correct segmentation in 83% of the cases. In 
the remainder of the cases, manual segmentation was required. Overall, these results 
5 are promising in that they offer a consistent, automated segmentation technique for the 
majority of juxtapieural nodules. 

Three-Dimensional Surface Generation - In addition to the three-dimensional 
voxel representation of a segmented nodule, it is useful to have a polygonal surface 
model for the segmented volume. Surface area and curvature calculations can be much 

10 more reliably performed on a smoothed surface representation. Computation of these 
metrics will be discussed in Chapter 6. 

Surface Tessellation - A closed 3D segmented region can be used to generate a 
polygonal surface model. One such surface model is a triangular tessellation of the 
surface. The surface models in this work are generated using a variant of the "marching 

15 cubes'" algorithm, described by Lorensen and Cline [48]. The modifications made to 
this algorithm will be explained, following a short description of the original 
implementation. 

In the standard marching cubes algorithm, each 2x2x2 voxel neighborhood is 
considered separately. The pattern of set voxels in this neighborhood falls in to one of 

20 15 basic classes. These are illustrated in Figure 35. Each of the 256 (2^) permutations ^ 
of set voxels in a 2 x 2 x 2 neighborhood can be generated using these classes via 
synmietry, rotation, and logical complement. 

For each of these canonical neighborhoods, a simple set of polygonal surfaces 
(composed of one or more triangles) is used to separate image background from 

25 foregroimd (see Figure 35). When these surface elements are generated for each 

neighborhood in a 3D image, a simple polygonal tessellation of the surface is created. 
The appropriate polygons corresponding to each 2x2x2 neighborhood may be 
assembled quickly using a lookup table indexed by a key corresponding to one of the 
256 possible neighborhoods, or octants, in the image. A description of the process is 

3 0 given in Algorithm 4.4. 
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Algorithm 4.4 (Three-Dimensional Surface Generation) 

for all iV^^,z = v(jc. . . jc + l^y.^.y-^ l,z.,.z + 1) 

compute octant o(x,y,z) 

index polygon look-up table using o(x,y,z) 
5 add triangles corresponding to o(x,y,z)j offset by ix,y,z) 

end 

Several modifications were made to the original marchiii^ cubes algorithm. 

First, the internal thresholding operation and gradient calculations have been removed. 

Segmentation decisions yielding the locations of each surface boimdary are made prior 
10 to polygonizadon, in the 3D segmentation method. The sxirface gmdient calculations 

needed for rendering Gouraud-shaded triangles are discarded, as the tessellation is iised 

for subsequent geometric analysis. Surface rendering is handled separately, by other 

tools in the VisionX system {\>rend, v3d). 

While the polygons described in the original algorithm are based on triangles, 
15 the effective aggregate polygons may be quadrilateral, or more complex non-planar 

polygons. In the new implementation, each of the constituent triangles is considered 

separately, so that the basic surface unit is the triangle, simplifying subsequent surface 

analysis. 

An additional, more interesting, change was made to the algorithm to avoid 
20 "holes" in the generated surface. When using the original marching cubes algorithm, it 
is possible for holes to appear in the final surface tessellation. It is possible to fix this 
problem by splitting several of the fifteen canonical classes and adding polygons to 
prevent these discontinuities [21]. 

One important note with regard to surface curvature estimation (described in 
25 Section 6) is that the triangles generated by the tessellation algorithm should be 

generated in a consistent right-or lefl-handed system. In other words, the order in which 
the three vertices of each triangle is encoded should be used to differentiate the inside 
fiom the outside of the surface. This will simplify surface normal calculations later in 
the curvature estimation method. 
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Surface Filtering - Although the surfece model of the segmented nodule is a 
much better approximation to the true nodule surface, it is still somewhat rough, as all 
angles are multiples of 45°. Algorithm 4.5 is a technique for jBltering the 3D surface 
representation. Each vertex in the surface is replaced with the weighted sum of 
5 neighboring vertices and itself. The parameter a specifies the proportional weight given 
to the neighboring vertices. When a = 0.9 , the new location of a vertex depends 90% 
on the locations of neighboring vertices and 10% on its original location. If a is set to 
1 .0, none of the original vertex location contributes to the result In addition, this 
smoothing algorithm may be iteratively applied n times, to control the degree of 
10 smoothing. 

Algorithm 4,5 (Three-Dimensional Surface Filtering) 

for ^= 1 : n 

for all Vf ^ V {for each vertex} 
Si 0 

15 for all e cidj(Vi) {for each vertex adjacent to Vi) 

Si{pc,y,z) ^ St + Vj 
end 

V,^ (1 - a) Vi + aiSil\adj{y;)\) 
end 

20 end 

The overall goal of the surface filtering process is to achieve a representation of 
the nodule surface that is more likely to represent its actual form. In combination with 
the isotropic resampling technique described in Section 3, it helps to mitigate the partial 
volume problem. One method of measuring the accuracy of the surface representations 
25 described is to measure the surface area of each representation of an object with known 
geometry. One such example is a sphere. 

A sphere of diameter 1 0 mm was generated synthetically and sampled to an 
isotropic grid with 1-mm resolution. Figure 36 shows shaded surface renderings of 
three surface representations of this sampled data. From right to left, they include: the 
30 origined voxel surfaces, the surface tessellation generated using the modified marching 
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cubes technique (see above), and the result of one iteration of Algorithm 4.5 applied to 

the previous representation. 



3 14. 15 mm^. The surface area calculated from these three representations are 480 mm , 

5 353, 16 nmi\ and 3 14.93 1 mm^, respectively. These correspond to surface area 
estimation errors of 52.8%, 12.4%, and 0.249%, respectively. Clearly, the filtered 
surface approximation gives the best estimate of the surface area. While the sphere may 
be an especially good case for the surface filtering algorithm, the surface of most 
nodules can be assumed to be locally spherical (in a small voxel neighborhood). In that 

10 model, the smoothed surface representation is a reasonably good one, and in any case, 
much more accurate than the voxel representation. Further discussion of surface area 
estimation is given in Section 6. The degree to which the surface is filtered, however, 
has an effect on all subsequent, surface-based metrics, including Ihe analysis of three- 
dimensional surface curvature distribution. This topic will be revisited in Section 6. 

1 5 Section 5: Volume Measurement -This section and much of the disclosure 

hereafter concerns volumetric measurements. The actual work and the discussion 
accompanying the work refers to "growth" in volumetric size. However, the present 
invention is not limited to "growth." It also includes the basis for measuring and 
handling information (e.g., data) related to reduction in volumetric size which might 

20 occur as, for example, by regression. Indeed, the scope of the invention covers "change 
and volumetric size. 

The most important step in volumetric measurement of ptilmonary nodules is 
segmentation of the nodule firom the surrounding lung parenchynaa and other structures. 
Segmentation methods were discussed in Section 4. Once the nodule region has been 

25 identified, volume measurement is straightforward. The volume in voxels, is merely the 
sums of all "set" voxels u in the region of interest (ROI). 



Analj^cally, the theoretical surface area for this sphere is equal to 1 OOtc mm 



M-l N-l Z-l 




(5.1) 



jc = 0 y = 0 2 = 0 
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To give the volume standard units, we must consider the volume of each voxel. 

M-1 N-l Z-1 
V(mm^ ) = S Z S z) • o,,; (5.2) 

:c = 0 y-0 z = 0 

In most CT data, the voxels are anisotropic, with much higher m-plane (x-y) than 
axial (z) resolution. The voxel volume is therefore: 

5 V„„i(m7W^)=Ux«ff •PywU«r = U^ •Vires (5.3) 

in - plane 

In our system, we resample the nodule ROI to isotropic space using triliaear 
interpolation. The voxel volume is then simply the cube of the isotropic resolution. , 

Ovoi(mm^)^vl^ (5.4) 

Isotropic resampUng, described in Section 3, has several benefits. It brings the 
10 data into isotropic space, where the resolution is the same in each dimension. This has 
the benefit of simplifying subsequent processing. Isotropic image filters may be used, 
rather than having to compensate for the different resolution in the less-resolved axial 
dimension. More important, the data are now at a supersampled resolution, where each 
original voxel value is responsible for several in the new image space. This helps 
1 5 mitigate partial volume effects in the original data. 

Exponential Growth Model - Exponential growth models have long been used 
in the differentiation of benign from malignant nodules [56, 89, 86, 80]. The basic model 
was developed on the assumption of uniform cell division, with cell dividing in two, in 
every generation. Therefore, the number of cells (N) as a function of the number of 
20 generations (n) is 

N(n) = Z 2' = = 2" --1 (5.5) 

This exponential growth can be characterized by a simple differential equation. 
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^ = XV (5.6) 
dt 



The solution of this eqmtion is straightforward. 



dt 



= XV (5.7) 



^ = Xdt (5.8) 
V 



]n.y=Xt+C (5.9) 



Solving for the constant C at r = 0; 

InFo = A-0 + C (5.10) 

C = InFo (5-11) 
Combining Equations 5.9 and 5.1 1 yields 

10 InV^Xt + lnVo (5.12) 

y = gCW+In Vo) = gdn TO) . (5_ J 3) 

This yields the most common form of the generalized exponential model for tumor 
volume based on an estimate of initial volume: 

V = Voe^ (5.14) 

1 5 where the ejqjonential coefficient X may be defined with respect to the nodule doubling 
time. 

X=hL (5.15) 
DT 
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A simple derivation follows, considering Vto be double Vo (this occurs at one doubling 
time). 

il = 2 = e^^ (5.16) 



0 



lii2 = ADT (5.17) 

A,=— (5.18) 
DT 

More generally, the exponential model can be expressed in terms of two volvime 
measurements At days apart: 

V2 



Vi 



= e^ (5.19) 



The doubling time, DT, of a nodule, then, can be computed using a ratio of two volume 
measurements and the time in days. At, between tiiem. This can be derived as follows: 



Fa = Vie*^ (5.20) 



— = e*^ (5.21) 
Vi 



ln(F2/Fi) = XdJ (5.22) 



Substituting our expression for X (Equation 5.15), results in 

ln(F2 / Vi) = ^"^'^ (5.23) 

DT 

m 

DT= ^^^-^ (5.24) 
ln(F2/Fi) 

Using Equation 5.24, doubling times can be computed for cliiucaliy-seen nodules, 

giv^ accurate measures of their volume at two times. Historically, however, the 
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doubling time has been estimated using measures of diameter on chest radiographs, or 
even in CT slices. In this case, there is an additional factor of three in the denominator 
(diameter varies as the cube root of volume). 

DT2D^ ^^^'^ (5.25) 
31n(D2/2)0 

5 To put the model in perspective. Figure 37 shows a comparison of characteris- 

tically malignant and benign doubling times as measured by the increase in nodule 
diameter from a baseline value of 1 mm. The time scale of the graph is 10 years. It can 
be seen that the progression of each nodule from the (currently) pre-detectable size of 1 
nam, varies significantly for the three doubling times, 75, 150, and 500 days. At the end 

10 of the 10 year period, the benign lesion is still only approximately 5 mm in diameter. 
The maUgnant lesion, however, reaches the T1-T2 size boimdary (30 mm) after only 6 
years. Furthermore, it can be seen that the malignant nodule remains in the small, sub- 
centimeter category for only 4 years ( » 1 500 days). It is worth noting that more 
aggressive malignancies, such as those with a doubling time of 75 days, will remain in 

15 the sub-centimeter category for only 2 years, and will reach 30 mm in diameter only one 
year later. 

Gompertzian Growth Model - Although the standard model for tumor growth is 
exponential, other models have been proposed. These models attempt to capture a 
possible retardation of growth rate for tumors of an advanced size, and have included 
20 logarithmic [25], logistic [57], and Gompertzian [45] mathematics as a means to model 
this effect. 

The latter model, based on a Gompertzian function, in which exponential growth 
is retarded at larger tumor sizes, is the most common of the altematives to the standard 
exponential growth model. This slowing of the growth rate has been hypothesized to 
25 correspond to cell death and central tumor necrosis, lack of nutritional supply, cell 
migration and differentiation, as well as homologous inhibition [1]. 

It is arguable how early in the history of lung cancer the exponential retardation is 
measurable, if present at all. It has been debated in the literature whether a Gompertzian 
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model is clinically practical [1], or observable in cancers outside those in animal models 
[89]. 

Much of the literature on lung cancer screening and growth analysis has been 
based on a simple exponential growth model. In addition, the focus of our study is small 
5 pulmonary nodules, which are likely at a stage prior to appreciable growth retardation; it 
has been reported that reductions from a constant growth rate may not be observable until 
the lesion is 2-3 cm in diameter [55]. These facts have lead to the choice of a simple 
exponential growth model for our characterization of nodule growth. 

InterScan Intervals - One of the issues surroimding doubling time estimation is 
1 0 the choice of InterScan interval; the number of days between CT examinations. The 

optimization of interscan interval, t\ is concerned with the following two constraints: (a) ti 
must be enough between scans so that accurate size change (grovs^) can be measured; 
and (b) ti must not be so long to allow substantial growth of a malignancy. 

Based on our experimental results regarding the accuracy of volume measurement 
15 as a function of nodule diameter d, (see Section 5), we can develop an expression for 
optimal values of the interscan interval. This leads to a scheduling function for early- 
repeat CT (ERCT) in the study of small pulmonary nodules. 

Our data on the reproducibility of volume measurements on synthetic nodules 
lead to the following assumptions. For small nodules (3 < d < 6), the percent error in 
20 volume estimation was 1.1% (RMS) and 2.8% (maximum). For larger nodules (6<d< 
1 1), the RMS and maximum percent errors were 0.5% and 0.9% respectively. If we take 
the maximum percent error, s for each size and conservatively estimate that the 
minimum reliably-detectable percent volume change, a is twice that (2 . e ), and, further, 
assume a conservative linear relationship between magnitude of error and nodule size 
25 (intuitively it should be cubic, decaying much faster with nodule diameter), we may 
construct the following relationship a as a function of nodule diameter, cf.Let s ^ and 

e^be the maximum percent errors in volume estimation for a small (d^ and large (di) 
diameter value, respectively. We assume a linear model, and that the minimum reliably- 
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detectable percent volume change, a, for a given diameter, d is 2Ed. This leads to the 
following model for a: 



a^md-^b (5.29) 



where the slope and intercept are 



2sL-2ss _ IjsL-es) 
di — ds di—ds 



(530) 



and 



b = 2es-mcis= ^^'^-^^ (5.31) 



di — ds 



The complete expression for a, then, is 



For a nodule of a given diameter, d, and doubling time, DT^ the minimimi 
reliably-detectable percent volume change is, a. given by Equation 532. We can then 
reformulate Equations 5.19 and 5.15 to yield the number of days needed to observe a 
percent change in volxmae of a, 

^ = l + (a/100) = ^^ (533) 

l + (a/100) = e^^2^^^^ (5,34) 

ln(l + (a /lOO)) = On2/DT)At (535) 

At = DT' ln(l + (a /100))/ln2 (5.36) 

This interval (At), can be used as the InterScan interval for nodviles of doubling 
time DT. It is the minimum duration in days needed to observe a reliably measurable 
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volume change of a percent, given a doubling time, DT. All that remains is the choice of 
doubling time on which to base our model. 

Doubling times for most malignant nodules range fix)m 30-400 days [47]. If we 
are conservative, we may estimate an upper bound on the doubling time of most 
5 malignancies as 500 days. We select a doubling time between that of the of the slowest 
growing nodule for which we would like to predict a status of meilignancy and that of the 
fastest growing nodule that we would like to classify as benign. If we tailor our InterScan 
interval to this doubling time, nodules with shorter doubling times (more aggressive) will 
be detected more easily, as they will have grown proportionally more in the time between 
10 CT examinations. Nodules with longer doubling times will exhibit little or no detectable 
growth in that time and will, therefore, be characterized as benign with respect to growth 
rate. 

Let DTo be the doubling time value that separates benign from malignant nodules. 
A conservative estimate of this value might be 500 days. We now have all of the 
15 mathematics in place to formally state the model for determination of the InterScan 



interval. 

Given a nodule seen at first HRCT to have a diameter of we may calcidate the 
value for the minimvim reliably-detectable percent volume change, a, using Equation 
5.32, and our estimates of and at diameters dg and di for the system under study. 

20 Following this calculation, the interscan interval (in days) for our chosen value of DTj^ is 
given by: 

t, = DT^ ln(l-h(a/100))/ln2 (5.37) 

Consider the following examples, based on the error estimates determined hereinabove, 
axxddiDTx) value of 500 days. From our synthetic nodule studies, « 3.0%, « 

25 1.0%, at ds = 3.0 mm and di = 11.0 mm. 
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EXAMPLE 1 

A 3.5-nim (diameter) nodule is detected using HRCT. Using Equation 5.32, we 
determine a to be 5 •75%. Equation 5.37 then yields a value of 40 days for tf. 

EXAMPLE 2 

5 A 9.0-mm (diameter) nodule is detected using HRCT. The corresponding value 

for a is 3.0%. The appropriate InterScan interval is 21 days. 

Each of these interscan intervals less than 90 days, perhaps the maximum time a 
second high-resolution CT study should be delayed. Intuitively, one must wait a longer 
interscan interval to accurately assess growth in a smaller nodule. As mentioned at the 
10 beginning of this section, one of the constraints of this model should be to avoid the 
possibility of significant nodule growth during the interscan interval. The maximimi 
amount of time a small pulmonary nodule (3 > d<10 mm) will be held between scans, 
given the error estimates described is approximately 42 days, which is less than 1.5 
doubling times for even the most aggressive malignancy. Traditional exam intervals for 
15 the study of lung cancer ranged from six months to two years (approximately 4-17 times 
as long), even though the smallest detectable lesion was at least 1 0 mm in diameter. 

This model for interscan interval calculation can be used as the scheduling 
function for duration to ERCT following an initial high-resolution study of a puhnonary 
nodule detected in a lung cancer screening program. While it may be administratively 
20 impossible to schedule patient visits on exactly the appropriate day suggested by the 

model, the value of ti serves as a lower boimd on the time needed to wait for an accurate 
growth assessment to be made. 

Accuracy of Doubling Time Estimates - An estimate the sensitivity of 

volimietric doubling time estimation can be made using e, the RMS percent error in 

25 volume measurement as determined experimentally (see Figure 38), For this model, as 
compared with the parameter a in the interscan interval model, it is sufficient to use the 
RMS percent error rather than the maximum error, as we are estimating the underlying 
variance of the measurements, rather than a conservative maximum error boimd, 
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15 



20 



The volumetric doubling time calculation is a non-linear function of two volumetric 
measurements which we will assume to be imcorrelated. In addition, it depends on the 
difference in time between scans. A/, which may introduce additional error. 



ln2 ■ At 



(5.38) 



The error in a nonlinear function, y —f (xj, X2, . . . , Xn), of uncorrelated variables can be 
estimated based on the variances of each variable. The exact differential of 



dx 



(5.39) 



n 



leads to the argument that the variance in the function may be estimated by a sum of the 
contributions of the variances of each variable. This is known as the law of error 
10 propagation. The variance in the nonlinear function above (if all jc/ are imcorrelated) is 



o\ + 



. dx . 



O ^ +... -I- 



^ 



K^nJ 



(5.40) 



In other words, the change (variance) in the function due to measurement error (as given 
by the exact differential) can be e>q)ressed as a function of the variance in measurement 
of each variable. We may, therefore, use this method to determine the sensitivity of the 
doubling time estimation, given estimates of the variance in each measurement, Fj, P2, 
and 

Starting &om Equation 5.38, we foim the exact differential of the DT function. 
The partial derivatives of jDT* with reject to each variable are 

SCDT) ln2 



a(A/) InCK^/F,) 



(5.41) 



dCpT)^ 

djPT) ^ -ln2 A/ 



(5.42) 



(5.43) 



and the exact differential of DT is therefore 



d(,DT) = 



Inl 



Inl' At 



V,(VniVJV,)) 



-1»2A/ 



F2(1«(F,/F,)) 



(5.44) 
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Then, using the law of error propagation, the variance in doubling tune estimation is 



ln2 



ln2 • tst 



( -In2-Ar 



(5.45) 



Interestingly, the variance of DT is not the appropriate error bound in this application. 
We would prefer, instead, to estimate the RMS error of the measurement, which is 
equivalent to the standard deviation, a, or the square root of the variance. 

For n measurements of variable x, these are equal to 



1 



RMS error = a = — 2 (x, - /i) 



N 



(5.46) 



1 = 0 



Thus, the RMS error in doubling time estimation is the square root qf Equation 5.45, or 



1«2 



\2 



i«(F2/n J 



1«2 • Lt 



N2 



V 



a\ + 



-In2-A/ 



V 



v^Mv^/Vi))^ 



(5.47) 



10 Doubling Time Estimation in Incidence Cases -Estimates of volumetric 

doubling time of a nodule can be made using two measurements of a nodule, as described 
by Equation 5.24. It may also be possible to estimate the doubling time of a nodule 
without two volxmie measurements, if there is a prior study in which the nodule was 
undetectable. Once a nodule is detected, prior CT examinations can be evaluated 

1 5 retrospectively to determine if a nodule had been present. If a nodule was in the pre- 
detectable phase on a prior CT, it may be considered an "incidence case," one that 
became detectable between regular examinations, as in a screening program. It is 
possible to estimate the maximum doubling time for these incidence cases using 
assumptions about minimum detectable size. 

20 If we define fi to be the minimxnn-detectable nod\ile diameter for a given CT 

protocol, it is possible to determine an upper bound on the nodule doubling time given 
two examinations; one in which the nodule was undetectable, and one in which precise 
volume measvirements were made. We may assume that the nodule size was smaller than 
P at the time of the first scan and compute the doubling time using the second study in 
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which the nodule volvime was measured and the interval between screening studies^ A/,, 
The actual doubling time of the nodule would have to be equal to or smaller (faster 
growth rate) than the derived value. The upper bound on the doubling time in an 
incidence case may be computed as 

^_Jn2^ (5.48) 

where Vfi'isa. minimiim-detectable volume estimate based on the minimum-detectable 
nodule diameter, fi. Given the assumption of a spherical nodule volume, would be 



(5.49) 



As an example, consider an aimual screening protocol in which 10-mm reconstructions 
10 are generated and for which the minimum-detectable nodule diameter is determined to be 
2.5 mm. 

EXAMPLE 3 

After several annual screens, a nodule is detected in a patient that was 
undetectable retrospectively in the previous exam. The nodule is scanned at high- 
1 5 resolution (1 mm slice thickness) for volumetric measurement yielding a volume estimate 
of 73,57 mm^. The upper bound on doubling time for this nodule could be estimated 
using = 8.18 mm^ and A^,, = 365.25. The incidence estimated doubling time would 
then be 

„ ln2- 365.25 . 
DT, = = 1 15.3 days 

' ln(73.57/8.18) 

20 Such a doubling time would put this nodule well within the range of malignancy. 

Furthermore, as this is an upper boimd on the doubling time, it is possible that the nodule 

is growing even more aggressively. 

Estimation of the doubling time of incidence cases may be an important 

altemative to waiting an appropriate interscan interval for a second volvimetric 
25 measurement. This is particularly true in incidence cases that are significantly larger than 

the minimum-detectable nodule size, as they are those most Ukely to be malignant. 
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EXAMPLE 4 

To validate the volumetric doubling time estimation method and to make 
comparisons with more traditional methods based on 2D measurement, a study was 
performed using high-resolution CT data acquired at two times for 13 patients enrolled in 
5 the ELCAP limg cancer screening study [28]. The final diagnostic status of their nodules 
was determined either by biopsy or demonstrated lack of growth for more than two years 
(2YNC). In addition, three nodules were included with unknown, but likely benign 
status. 

The 3D resampling and segmentation procedures described in this work were 
1 0 used prior to automated measurement of nodule size. Volumetric measurements were 
made of each scan of each nodule. For comparison, area meeisurements were made using 
the CT slice of maximimi cross-sectional area, and diameter esthnates computed from the 
major and minor principal axes of the nodule in the same image. 

Nodule doubling times were computed for each case based on two criteria; 
15 volumetric measurements (Equation 5.24) and area measurements (Equation 5.25). Table 
5.3 Usts the diameter, volume, and area estimates for each case, as well as the two 
doubling time determinations and the final diagnostic status. 

Volumetrically-determined doubling times for the malignant nodules ranged from . 
51 to 177 days, all well within the range for malignancy. Similarly, doubling times for 
20 the non-malignant nodules ranged from 3 1 8 to 33700 days for growing nodules and from 
-242 to -3610 for those nodules that were found to have reduced volume on repeat 
examination. With the exception of Case 14, all non-malignant nodules had doubling 
times consistent with benignity. 

An example illustrating volunietric growth estimation is shown in Figures 39 and 
25 40. Figure 39 shows 2D CT images of a small pvdmonary nodule (baseline diameter z: 
5.5 mm) taken at baseline and on repeat HRCT 33 days later. The marked growth of this 
small nodule is clearly visible even with such a short InterScan period. The nodule data 
were resampled to 0.25 mm isotropic space and segmented using Algorithm 4.2, which 
performs morphological segmentation to remove the small vascular attachments. Figure 
30 40 shows surface-rendered 3D representations of the segmented nodule at baseline and on 
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repeat HRCT, respectively. The overall growth of this nodule is even more evident in the 
3D representations than from the original 2D CT slices. Note the asjrmmetric growth 
observable in the exaggerated protrusions on the top and front nodule surfaces. 

Volumetric analysis was performed on both scans of the nodule. It was foxmd that 
the nodule grew from 62.5 mm^ to 85.3 mm^ during the InterScan period. Using the 
exponential nodule growth model and Equation 5.24, the nodule doubling time 
Table 5.3: In-vivo nodule doubling times based on change in volume and area 
measurements 
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was computed to be 74 days, well within the range for malignancy. The nodule was 
subsequentiy resected and proven malignant by histologic examination. 
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The asymmetry observed in this example was relatively minor. Non-miiform 
growth in three dimensions may, however, be a significant factor in doubling time 
estimation. When comparing volumetrically-determined doubling times with those based 
on 2D area measurements, we see discrepancies is several cases. Consider Case 1, where 
5 the volumetrically-determined doubling time is markedly different firom that based on 
area measurement Figure 41 shows 2D images of the CT images vwth maximum cross- 
sectional area used for computation of area-based doubling time. 

The segmented images of the nodule at baseline and at repeat HRCT are quite 
similar. Their cross-sectional areas are 36.5 and 36.6 mm , respectively, and their 

10 perimeters are 22.7 and 23.4 mm, respectively. From this two-dimensional perspective, it 
is impossible to discern that the nodule is growing in any sigmficant manner. Figure 42, 
however, illustrates tiie nodule in 3D, viewed perpendicular to the scanner axis. This 
reveals that the nodule is, in fact, growing perpendicular to the scan plane. Thus, while 
the area-determined doubling time for this nodxile was a benign-appearing 9700 days, the 

1 5 true volumetrically-determined doubling time was computed to be 104 days, consistent 
with malignancy. Again, the nodule was resected and confirmed malignant by histology. 
Similar cases in which 2D measurements were insufficient to correctly characterize 
growth were Cases 7, 9, and 1 1 . True volumetrically-determined doubling times appear, 
therefore, to be better quantifiers of nodule growth, as would be expected. 

20 Volumetric Growth Index - When evaluating the statistical differences between 

doubling times computed for malignant and non-maUgnant nodules, it becomes apparent 
that the noiilinearity of the measure poses a challenge to evaluation. Students t-test was 
used to estimate the difference in means between the V£dues of DT computed for 
malignant and non-malignant nodules. The fact that both large doubling times 

25 (corresponding to slow growth) and negative doubling times (corresponding to reduction 
in nodule size) both correspond to non-malignant processes, assessment of tiie means of 
the DT distributions does not effectively characterize their differences. 

To eliminate this non-linearity and to provide another intuitive measure of 
growth, we may define the volumetric growth index, or VGL The VGI is defined to be 

30 the relative increase in nodule volxmie in one year. This effectively remaps the non- 
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malignant nodules that grow slowly or that are seen to be reduced in volume to the range 
of small positive and negative numbers. This transformation can be accomplished as 
follows. Beginning with Eqxiation 5.24, substituting one year (365.25 days) as Ar, and 
solving for the ratio of volume measurements we get 



— = exp 



rin2- 365.25 ^ 



DT 

Subtracting one from this ratio, to give the relative increase in volume yields the 
expression for VGI: 



(5.50) 



VGI 



<ln2- 365.25 ^ , /c ct\ 
—1 C5.51) 
DT ) 



This metric describes the fractional relative increase in nodule volume each year, given in 
10 units of relative volume (percent/100). Consider the following three examples. 

EXAMPLES 

A benign nodule is determined to have a volumetric doubling time, DT, of 1240 
days. Using Equation 5.51, we determine the volumetric growth index, FG/, to be 0.227, 
or 22.7%. This indicates that at the cvuirent growth rate, the nodvde will increase 22.7% in 
15 volume per year. 

EXAMPLE 6 

A malignant nodule is determined to have DT value of 75 days. We may then 
determine its volumetric growth index, KG/, to be 28.2, or 2820%. In each year, this 
aggressive nodule with increase its volume approximately 28-fold. 
20 EXAMPLE 7 

A benign nodule is determined to have DT value of -21 00 days. The 
corresponding value for KG/ would be -0.1 14, or -1 1.4%. The nodule is therefore slowly 
decreasing in size. 

EXAMPLES 

25 With the volimietric growth index, we may now perform a statistical comparison 

of the distribution of this metric in maUgnant and non-malignant nodules. Table 5.4 lists 
the doubling time and volumetric growth index for each of the in-vivo nod\iles given in 
Table 5.3. 
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Statistical analysis reveals that the distributions of KG/ values for malignant and 
non-malignant nodules are significantly different The VGI for malignant nodules was 
40.44 ± 25.53 (mean ± sem), whereas for non-malignant nodules it was 0.2620 ± 0.1520. 
The distributions showed a significant difference in means, using Student's t-test (p < 
5 0.028), Given the limited sample size, and to remove assumptions about the normality of 
the distributions, we may also use non-parametric statistics to compare values for the two 
classes of nodules. Using Wilcoxon rank sums, the distributions were also shown to be 
significantly different (p < 0.002). 

Table 5.4: In- vivo volumetric doubling times and associated volumetric growth indices 
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1 0 Section 6: Analysis of Three-Dimensional Moments - Three-dimensional 

moments [65] have been used to characterize shape information in a variety of three- 



* 
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dimensional data. The standard 3D Cartesian moment set of order m(pqr), contains all 
moments m(pqr), such that p+q-^r < n. These moments are defined as follows: 

CO CO CO 

pgr =^ 



= J J J x''y^z''f(Xyy,z)dx dy dz 



— 00 —00 



(6.1) 



where / (jc, y^ z) is a continuous function of three dimensions. In 3D sampled images 
(such as those arising in from CT data) this results in a discrete image v (x, y^ z) of size 
(MxNxL). The corresponding moment set can then be calculated using the discrete 
expression 

= S S S x''y'z'vix,y,z) (6.2) 

x=0 >=0 2=0 



rn 



Many of the lower order moments have physical interpretations. The zeroth- 
10 order moment in three dimensions, for example, is the volume of the object: 

M-l N-l L-l 



^"000 =2 S S v(x,j;,z) (6.3) 



;c=0 ;;=0 2=0 

Similarly, the center of mass of the object is described by (x, y^ z), where 



x^::^y^:::^z=:=^^ (6.4) 

^000 '^000 ^000 

Geometric and Densitometric Moments - Moment analysis methods for the 

15 characterization of objects in images fall into two general classes, based on tjrpe of 

application: applications that are concemed purely with shape characteristics, and those 

that take into consideration the distribution of intensity within regions. These lead to 

two types of moment sets that can be computed as follows: (a) geometric moments - 

each voxel value is considered as binary-valued, contributing a 0 or 1 to the moment 

20 sum; and (b) densitometric moments - each voxel value is considered to have a range 

of intensities, or densities. 

Geometric moments are computed from a segmented, binary representation of an 

object or scene, whereas densitometric moments are computed from a gray-level 

representation. Consider Equation 6.2. When computing geometric moments, each 

25 value of v (x, y, z) is binary-valued. In contrast, when computing densitometric 
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moments, v(x, yy z) may have a range of values in a suitable representation (8-bit/ 16-bit 
integer, 16-bit/32-bit floatingpoint). 

Both geometric and densitometric moment analyses can be used in the 
characterization of pulmonary nodules. These techniques will be discussed in Sections 
5 set forth below. 

Moment Invariants - In an effort to develop metrics that do not change with 
object scale, position, or orientation, certain nonlinear combinations of moments have 
been devised known as moment invariants. These are of obvious use in both pattem 
recognition and machine vision, as they allow objects to be described in a unique 
10 manner. One drawback of conventional moment invariants is that they have limited 
physical meaning (as do the higher-order moments themselves). 

Invariant Moment Sets - An alternative technique in object description is to 
establish a standard set of moments that can uniquely describe an object. It can be 
shown that such a unique set would need to be of infinite order, but that moment sets of 
1 5 order 3-5 are sufficient for practical classification of most three-dimensional objects 
[65, 60]. 

The first step toward generating an invariant moment set is to require the center 
of mass to be located at the origin. 

j"ioo = /^oio - /^ooi = 0 (6.5) 
20 This produces a central moment set, {//pgr}- If further normalizing restrictions 

are placed on a moment set {mpgr}y it can be transformed into a standard moment set, 
{Mpgr }. Briefly, these restrictions are as follows: 

1. Volume-The object volume must be normalized to 1 .0. 

25 Afooo=1.0 (6.6) 

2. Center ofMass-Thsi center of mass of the object must be at the origin. 

Mm = Mow = Mqoi = 0 (6. 7) 

3. Principal Axes-Tho^ principal axes of the object's ellipsoid of inertia (EOI) 

must be aligned with the Cartesian axes. 

30 Mno = Mioi = Mon = 0 (6.8) 
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4. Moments Inertia-The moments of inertia must be (in descending order) z, y. 



andx. 



(6.9) 



5. Rotational Normalization-ThQrc must be a unique orientation of the object 



W300 ^ 0, Wo3o ^ 0, and (optionally) m^^^ > 0 



(6.10) 



Transformations of Three-Dimensional Moments -Consider the following 3- 
by-3 matrix, specifying a generalized orthogonal transformation of the 3D image space. 



'^Moo Moi /^02 ^ 
A^io A^ii f^i2 

^^20 M2I ^^22 J 



(6,11) 



10 



15 



20 



The efifect of this arbitrary orthogonal image transformation is to transform the 
moment mpgr to /Wp^r- Its effect can be achieved directiy in moment space using the 
following formula: 



I7t = 



p q r 

E E ■ Z Z E 1: C U m,,, (6.12) 

5j=0 =0 ^2 =0^2=0 ^^3=0 ^3=0 



where the combinatorial product C is defined as 



P 



V2 J\^3j 



the transformation product U is defined as 

and the indices a, by and c of the independent moment term are as follows. 

a = ^1 H- /2 + ^3 
b = Si+S2'^S3-tj''t2't3 



(6.13) 



(6.14) 



(6.15) 
(6.16) 
(6.17) 



This formulation allows scaling and rotation to be applied in one step. A 

complete a^uie transformation (translation, rotation, and scaling) and may be 

86 



wo 01/78005 



PCT/USOl/11820 



accomplished by translation of the center of mass to the origin followed by a linear 
transformation combining scaling and rotation. 

Object Orientation -In order to put the object under study into standard form (and 
thereby compute {Mpqr})> it is necessary to determine its orientation in three-dimensional 
space. This can be accomplished via the solution to the following eigenproblem: 

Ax=-%x (6.18) 



where 



^no ^020 ^on 



(6.19) 

The orthononnal basis of eigenvectors produced in the solution of this problem 
1 0 will point in the directions of each of the principal axes of the object Using these unit 
vectors K^, Ky, Fz, it is possible to find the orientation of the object. Consider a standard 
orientation where the major principal axis (f^,) is aligned with the x-axis, the 
intermediate principal axis (V^ is aligned with the y-axis, and the minor principal axis 
(F^) is aligned with the z-axis (Figure 43). 

1 5 Note that the three principal axes must be mutually orthogonal, as they form an 

orthononnal basis for the object distribution. This standard orientation can also be 
specified in terms of second order moments (moments of inertia): 

nttm ^ ^^020 ^ ^002 (6.20) 
The moment of mertia about the x-axis (maoo) should be greater than or equal to 
20 that about the y-axis (wo2o), which should, in turn, be greater than that about the x-axis 

(/W200)- 

The lengths of the principal axes (twice the length of the semi-axes of the 
ellipsoid of inertia) can be derived fi:om the eigenvalues of the system. If the eigenvalues 
of the system given in Equation 6. 1 8 are sorted such that 
25 A.o>^i> (6.21) 

then the lengths of the principal axes are 
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length=K| = 2VI7- j f 



(6.22) 



(6.23) 



height = \V^\ = 2ylX^ 



1 3K 



(6.24) 



Given the restrictions of standard orientation, it is possible to define notions of 
roil, pitch, and yaw, allowing the specification of an arbitrary orientation as three 
rotations firom the standard form. If we restrict objects to be longest along their major 
principal axis and shortest along their minor axis, these definitions correspond to their 
standard geometric counterparts. Simply, roll is a rotation about the major principal 
axis; pitch is a rotation about the intermediate principal axis; yaw is a rotation about the 
minor principal axis. Using the principal axis vectors obtained firom the solution to the 
above eigenproblem, Ihese rotational angles are simple functions of their projections on 
the KZ, AZ, and A7 planes, respectively. 



roll = cos 



-1 



pitch = cos 



V 



-1 



sign(F^ (y)) 



(6.25) 



J 



\ 



sign(K,(z)) 



(6.26) 



yaw = cos 



-1 



•sign(p;();)) 



(6.27) 



The projections, derived from the dot products of each vector with the line of 
minimum distance in each corresponding plane simplify to the following: 

projy^V, = (0. V^), F,)z)) (6.28) 
proj^V:, = (F:c(x),0,F,)z)) (6.29) 
proJ^V^ = (Fx(jc),FxO'),0) (6.30) 
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The sign 0 (signiim) function is used to uniquely define the rotations, as there are 
otherwise two possible solutions mapped to the same angle by the inverse cosine (the 
range of cos'^ is [0,7i;]). 

Table 6.1 : Primary size metrics in two and three dimensions 



Feature. 


2D 


3D 


gross size 


area 


volume 


surface 


perimeter 


surface area 


dimensions of EOI 


length, width 


length, width, height 



5 

Once the orientation of an object has been determined, it is possible to 
transform the entire moment set {jtipgr} into its corresponding standard moment set, 
{mpqr}p using the general transformation given in Equation 6.12. 

Primary Metrics -Primary metrics of nodule size and density are those that are 
10 computed directiy from the segmented nodule data. 

Primary Size Metrics - The primary measures of nodule size in two and three 
dimensions are shown in Table 6.1, where "EOI" refers to the ellipse of inertia (in 2D) or 
ellipsoid of inertia (in 3D) representing the nodxile. Calculation of two-dunensional 
metrics is normally performed on the single CT image of maximum cross-sectional area. 
15 Volume and area calculations are straightforward, as they are the zeroth-order 

geometric moment in 3D, and 2D, respectively. Using the general equation for 
moments (Equation 6.2), we get volume (/wooo) as expressed in Equation 6.3. Similarly, 
the area of a two-dimensional object is moo* It is important to note that these equations 
result in size meeisures in units of voxels or pixels. Section 5 discusses the 
20 normalization of volume imits based on voxel resolution. 

Surface Area Calculation - There are two types of methods for surface area 
calculation. Those that operate on segmented voxel data, and those that operate on a 
polygonal surface representation of the nodule (described in Section 4 Surface area 
calculations made directiy from the segmented voxel representation are basically the sum 
25 of the exposed voxel faces, those that separate object from background. The area of each 

exposed voxel face is simply the product of two of the image resolution values. 
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Algorithm 6.1 (Surface Area Calculation (Voxel-Based)) 

for all v{x,y^z) b (y{x^y^z) = l) {for each foreground voxel} • 
if (v(jc - 1, z) = 0||v(x + 1, y, z) = o) 
5 <- 5 + (v^ • v^) {add yz rectangle} 

if (v(x, y-U)=^ 0||v(x, >; + 1, z) = O) 

5 <- 5 + (v^ • v^) {add xz rectangle} 
if (v(x, J/, z - 1) = 0||v(x, j^, z + 1) = O) 
S <-S + {y^'Vy^) {add rectangle} 

10 end 

This voxel-based algorithm tends to exaggerate the surface area estimate, as the 
rectangular approximation to the surface is, in general, not nearly as smooth as the 
object it represents. Figures 44 through 47 illustrate this point in two dimensions. 

Figure 44 shows the silhouette of an object whose shape will be discretized by 

1 5 the sampling procedure. Figure 45 shows the discretized representation of the object. 
Note the error introduced by sampling the continuous image to the grid. Similar to the 
partial volume problem in CT, pixels containing insufficient overlap with the object are 
assigned to the backgroimd. The boundary of this pixel-based representation clearly 
overestimates the true extent of the object surface (perimeter). Figure 46 illustrates a 

20 polyline representation of the object boundary. The polyline (a two-dimensional analog 
to the polygonal surfaces used in our nodule analysis) provides a better representation of 
the true object form and corresponding surface length. Finally, Figure 47 illustrates a 
filtered (smoothed) polyline representation, where the location of each vertex is replaced 
by a weighted sum of its original location and that of the two adjacent vertices. 

25 Figure 48 illustrates the difference in voxel and surface representations in three 

dimensions. Each is a shaded surface rendering of a polygonal surface representation 
of the same small (d - 5 mm) pidmonary nodule. The image on the left is the surface 
given by the voxels themselves. The center image is the result of surface tessellation by 
the modified marching cubes algorithm discussed in Section 4. The image on the right 

30 is the filtered version of the center image following two iterations of Algorithm 4.6. 
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It is clear from the figure that the voxel-based representation exaggerates the 
surface area. The degree to which the surface area is overestimated can be evaluated by 
comparing voxel-based surface area estimates to those computed from polygonally 
tessellated surfaces. 

The standard sxu&ce area estimation method for this work operates on the 
polygonal surface representations of the nodxile described in Section 4. Surface area is 
computed as the sum of areas of each polygon in the svuface. Algorithm 6,2 is a 
description of this method for surfaces tessellated with triangles (the standard). 
Extension of the method to include convex planar polygons of n sides is straight- 
forward, as they can be decomposed into triangles (as is the case for the squares used to 
tessellated the voxel-based image). It is based on computing the area of each triangle 
using Heron's formula. For clarity, we define the Euclidean distance function between 
vertices Va and Vb as: 

D{y.>v,)=4{vXx)-v,{x))U{vM-vMY +{v,{z)-v,{z))' (6.31) 

Algorithm 6.2 (Surface Area Calculation (Polygon-Based)) 

for all Ti= {VcVuVa} {for each triangle} 
a = D{V„V,) 

s = —(a + b + c\ 

S<-S + A 
end 

Consider the nodule surfaces shown in Figure 48. Using Algorithm 6.2, we find 
the surface areas are 122.25 mm^ 81.2759 mm^ and 79.4623 nml^ for the voxel-based, 
polygonal surface, and filtered polygonal surface representations, respectively. 

This result is similar to that derived in Section 3, where it was shown that the 
error in area and volume measurement is reduced with an increasing resampling ratio. 
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Filtering of the polygonal representation provides much the same benefit, allowing the 
locations of each vertex in the boundary to have displacements that are more finely 
resolved than on the original grid. Further evidence of the benefit of surface filtering 
prior to the estimation of surface area was shown in Section 4. The error in surface 
area estimation of a 10 mm sphere was shown to be reduced firom 12% when calculated 
from a polygonal surface generated by the modified marching cubes method to 0.25% 
following the application of Algorithm 4.5. 

An additional experiment was performed to assess the relationship between the 
nimiber of iterations of surface filtering and estimation of surface area. Synthetic 
spheres of varying diameter (from 2.5 to 25 mm) were generated, tessellated, and 
surface filtered a varjdng number of times prior to surface area computation. The 
surface area estimates were within 2% (mean = 0.4%) of theoretical values following 
one filter iteration as compared with 12% (mean = 9.4%) using the imfiltered 
tessellations. Furthermore, for spheres greater than 5mm in diameter, the surface area 
estimates could be improved to within 1% (means == 0.4%) following two iterations of 
Algorithm 4.5. 

Primary Density Metrics - Computation of voxel density statistics are based on 
statisticEd moments. More specifically, they are functions of the central statistical 
moments. Analogous to the central moments of the image data described in Section 6.2, 
these moments are summations of powers of the voxel density values, normalized to the 
mean value. 



This is exactiy analogous to area and volume in two and three dimensions, 
respectively. Note that the nomialization of these moments is to the mean value in the 
one-dimensional distribution of density, not to a particular spatial location (e.g. the 
center of mass). Therefore, voxel density statistics provide information about the 



0 



(6.32) 




(6.33) 
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Statistical distribution of all voxel densities in the nodule region without regard to 
spatial location. 

The mean voxel density is the arithmetic mean of all density values in the 
segmented nodule volimie and can be defined as the first-order moment divided by 
5 zeroth-order moment 

fi=— S o{x,y,z) (6.34) 
0 

Variance, a measure of the variability of density values, is derived firom the 
second-order moment, expressing the mean-squared deviation of the voxel values firom 
the mean. 

10 variance = — = — z](v(x,>',2) — /x)^ (6.35) 

An xmbiased estimate of variance is sometimes used when it is hypothesized that 
the sample mean does not represent the true mean of a variable over an entire 
population. Such hypothesis are not usually necessary when dealing with the 
distribution of voxel values in a nodule, as we often have data representing the entire 
1 5 nodule volume. 

variance(unbiased) = — — = XK^C^* l^Y (6.36) 



The standard deviation, a is another measure of the variation in values firom the 
mean. It is simply defined as the square root of the variance. 

<j = V variance (6.37) 
20 Two additional parameters based on statistical moments, skewness and kurtosis, 

are used to quantify the shape of the distribution. Skewness measures the shift of the 
density distribution to the right or lefi: (above or below) the mean. 

Kurtosis is a measure of the "peskiness" of the distribution. 
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kurtosis - ^ - 3 = ^ ^ ^ ^ - 3 (6.39) 

Note that the expression given above is normalized (by subtracting 3 from the "raw" 
kurtosis) such that the kurtosis of a normal distribution is 0.0. 

These voxel statistics can be used to assess the uniformity of density in the 
5 nodule, which may be useful in characterizing nodule cavitation, or, more generally, 
cell necrosis. In addition, the mean voxel density within the nodule can be used in the 
assessment of nodule calcification, a major predictor of benignity. If the mean density 
or a significant subset of the density distribution is found to be higher than that 
observed for calcium, the likelihood of benignity is considerably higher. Patterns of 

10 calcification within nodules have been studied [47], but are generally only observable 
in nodules larger than 1 cm in diameter. In this study, however, our focus is generally 
restricted to small non-calcified nodules. 

Secondary Metrics - Secondary metrics of nodule size, shape, and density are 
those that are derived from primary metrics, rather than meeisured directly from the 

1 5 image data. One of the major benefits of secondary metrics is that they are typically 
defined to be size invariant Unlike primary size metrics, the size invariance of 
secondary metrics allows nodules of different sizes to be compared without regard to 
size. Although overall nodule size is a metric typically predictive of malignancy, it is 
less appropriate as a measure of small pulmonary nodules. Size as a predictor of 

20 malignancy is typically used to stratify larger lung masses lung masses (d>3 cm) (not 
generally classified as nodules) firom smaller lesions. These large masses have a high 
probability of malignancy, as benign nodules rarely grow to this size (82). In the study 
of small pulmonary nodules (^ < 1 cm), however, such differentiation based on size is 
not possible. Therefore, the study of size-invariant metrics is concerned with 

25 detenniiiing those shape and density characteristics that can be used to differentiate 

benign firom malignant lesions when they are all of a similar size. 

Aspect Ratios - The simplest secondary metrics describing the shape of a nodule 

are the aspect ratios, simple ratios of the dimensions of the segmented nodule volume. 

These dimensions are computed using moment analysis and are (in descending order of 
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magnitude) the length, width, and height of the ellipsoid of inertia describing the 
segmented volume. Computation of these dimensions is done using Equations 6.22 - 
6.24, following the solution of the eigenproblem described above. The three aspect ratios 
of these nodule dimensions are defined as 

iH8 = f^ (6.40) 

height 



LWR = ^ (6.41) 
width 



WHR^^^^!^ (6.42) 
height 

10 Similarly, an aspect ratio may be defined in two-dimensions based on the ellipse 

of inertia computed from a single 2D image. In this case there is only a single ratio, 
that of length (major axis) to width (minor axis). 

Nodule Compactness - Another important measure of the geometric distribution of 
a pulmonary nodule is compactness. Compactness is a measure of the ratio of size to 
1 5 surface of a nodule. 

In two dimensions, compactness has been defined as a function of cross- 
sectional area and perimeter (22). 

^ ^ 4n'Area 

Compactness2D = (6.43) 

Perimeter 

The constant, Aln, and exponent of perimeter are used as a normalization such 
20 that the 2D compactness of a circle is equal to one. 

We can extend compactness to a three-dimensional metric. Three-dimensional 
compactness may be defined as a fijnction of nodule volume and surface area [62]. 

^ 6yf7t'V 

Compactness = (6 .44) 

S 

The constant, S-y/TV , and exponent of surface area are used to normalize the 
25 fimction such that the compactness of a sphere is equal to one. 
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Sphericity and Circularity - Two additional measures of nodule shape may be 
defined that capture both the compactness and the major;minor aspect ratio the nodule. 
Sphericity is a measure of similarity of the nodule to a sphere of the same volume. 

Spericity^^^^^^^^^^^^ (6.45) 

LHR 

5 In two dimensions, the analogous metric is circularity. 

Circularity ^ ^'''^P'''''^^''^- (6.46) 

LWR 

Secondary Density Metrics - Secondary metrics of density distribution within 
the nodule may also be defined. The goal of these measures is to quantify the regularity 
of the nodule density distribution and how it deviates firom a imiformly dense sphere. 
10 Eccentricity, ^5 of the density distribution measures the displacement between 

the geometric and densitometric centers of mass (COM). It is defined as 

B, = d(COM^COM^) (6.47). 

using the D operator for Euclidean distance (Equation 6.3 1^. To obtain a truly size- 
invariant metric, the eccentricity may be norm2ilized by the cube root of the volimie (an 
1 5 estimate of nodxile radius). 

I ijCOM^OM^) ^^^^^ 

Density skew, 0^ is a measure of the angle between the geometric and 
densitometric ellipsoids of inertia (EOI). If we define the orientation of the geometric 
EOI to be dgeom and that of the densitometric EOI to be Ociens> the density skew is simply 

20 <l>,=e^.-e^. (6.49) 

One concern in using this metric is the stability of the measure of EOI 
orientation. Note that as the nodule becomes more spherical, the calculation of the 
orientation becomes xmstable. Therefore, the measure of density skew is less accurate 
as the nodule is more spherical. This leads to the use of the normalized density skew, 
25 0ci » which is simply resealed by the sphericity of the nodule such that the less- 
confident values are closer to zero and the more-confident values are given a greater 
magnitude. 
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e -e. 

A ^ geom, dens. 

^ Sphericity ^ ' ^ 

Curvature Model - Analysis of surface curvature is the study of the rate of 
change of the surface normal, 4>, with respect to the surface length- In two dimensions, 
5 this is the derivative of the normal vector with respect to the arc length. 

k:=-^ (6.51) 
as 

In three dimensions, an infinite niunber of curvatures can be defined at a point 
on a surface, as they are computed from the curve given by the intersection of the 3D 
surface with any plane containing that point and that is parallel to the surface normal at 
1 0 that point 

In our discrete piecewise linear model of the 3D surface of a nodule, the 
situation is more straightforward. The surface curvature can be computed as the 
change in surface normal between a particular vertex and any (or all) of the adjacent 
vertices. Curvature estimates may be defined at each vertex, or on each triangle in the 
1 5 surface. Methods for the computation of these curvature estimates will now be 
described. 

Curvature Estimation Method - Curvature estimation methods are often based 
on analysis of 2D and 3D image gradient information. In this work, curvature 
estimation is performed on a smoothed piecewise linear surface model, derived fi-om the 
20 results of 3D image segmentation. The surface tessellation and smoothing are described 
in Section 4. 

Normal and Curvature Estimation — Three-dimensional surface curvature 
estimation may be perfomied on the nodule sxjrface model in two steps. Curvature 
estimates are made for each vertex comprising the surface and, in turn, for each 
25 triangle. 

For each vertex Vi in F, connectivity analysis is performed to generate an 
adjacency Ust, adj(V^. From this adjacency list, it is possible to compute the set of 
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triangles, T, of which Vt is a member. Consider a pair of vertices (K^, Fb) adjacent to Fi. 
If they also share adjacency, there exists a triangle with vertices {Fi, F^, F^}. 
This method is described more foraially in Algorithm 6.3. 

Algorithm 6,3 (Triangle Detection) 

VF, eF 

V(F^,F^) e adjXVf) {For each pair of vertices adjacent to F,} 
if F, e adjXV,) {If a4r(FJ contains FJ 
r = (ru{F„F„FJ) 

end 

Figure 49 illustrates a small patch of a 3D tessellated surface. Vertex F/, and 
the four vertices in adj (V^ = {Va, F& Vc V4} are shoA?vn. This patch will form the basis 
of a running example to be used in illustrating the curvature estimation model. 

Once the set of triangles, r= {To . . - Tm}y of which Fj is a member has been 
determined, the first step is to calculate the surface normal for each triangle in T. If T2 
= {Fi, Va, Vb}\ then its surface normal,^^i, can be computed as the normalized cross 
product of two of its sides: 



(6.52) 



This relationship is illustrated in Figure 50. Note that we normalize this result to 
obtain the xmit surface normal vector for each triangle. The surface normal at vertex Fi, 
0u can now be calculated as the average of surface normals of each triangle of which it 
is a member. 

m 

<l>^-^(Z^j)^\T\ (6-53) 

This is illustrated in Figure 5 1 for the vertex V. The curvature at that vertex, Fi, is 
computed using the four triangles of which it is a member, {Ta, T^, Tcy 7^}, and their 
respective surface normals, {y/a y/b Wc y^d}- 

Curvature estimates for each vertex, C, may be derived from the vertex surface 
normals, <t>. The goal is to estimate the derivative of the sxirface normal across points 
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on the surface. Since the tessellated surface model is a piecewise linear approximation 
to this surface, discrete curvature estimates may be made at each vertex by considering 
the variation in vertex surface nomial and that of adjacent vertices. Such an 
arrangement is illustrated in Figure 52, which shows the surface normals at Vi (4> /) and 

* 

5 each of the vertices in adj (Vi). Note that the surface normals of the adjacent surfaces 
have been determined using triangles not shown in the figure. 

Based on this foxmdation of triangle and vertex surface normals, several types of 
curvature estimates may be defined. Each of these is based on the angular difference 
between the surface normal of a vertex, F/, and that of each of the adjacent vertices, 
10 adjXVi). The angular difference, ft,, between surface normals at vertices Fi and Va can be 
computed iising the arc cosine of the normalized scalar (dot) product of the two surface 
normal vectors, and 



<l>„ = cos' 



(6.54) 



Graphically, each of the vertices (and associated surface normals) may be translated so 
15 that they are coincident, and the angle between there clearly defined. This is illustrated 
in Figure 53. 

With this differencing mechanism in place, we may now discuss the varying 
curvature estimates possible, given our 3D tessellated model. The simplest of these is 
the determination of the principal curvatures. The principal ciurvatures, ki and K2 are 
20 defined as the maximum and minimum values, respectively, of the normal curvature 

evaluated at a vertex. Note that in our discrete vertex space, this may be approximated 
using the extrema of the curvature values with respect to each of the adjacent vertices. 
Algorithm 6.4 illustrates the computation of principal curvatures for each vertex in a 
3D surface model. 

25 
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Algorithm 6.4 (Principal Curvature Estimation) 

for all e V{for each vertex) 
for aHV^^adJiV,) 

0, =cos J 
end 

fc^(F;.)=max(0) 
K:2(F,)=min(0) 
end 

Two additional curvature measures may be defined in terms of the principal 
curvatures, ki and K2. Gaussian cvirvature, is the product of the principal curvatures 
5 at a vertex. 

K^K^K^ (6.55) 
Mean curvature, H, is the mean of the two principal curvatures. 

Gaussian and mean curvatures may be computed as a simple extension to Algorithm 
10 6.4. While mean curvature is an excellent metric for continuous surfaces in 3D space, 
it is much more sensitive to noise in oxir discrete model. Therefore, a more robust 
measure of the average normal curvature at each vertex is used. Rather than computing 
if based on only two atomic curvature values, we may estimate the average normal 
curvature; C, as the average of angular differences in surface normal between a vertex 
15 and all of the adjacent vertices. This method is described in Algorithm J. 

Algorithm 6,5 (Average Normal Curvature Estimation) 

for all Vf (for each vertex) 

(compute curvature at vertex as average of angular differences) 
for all V € adj'iV.) 

0,=cos-'((^,.0J/(j</»,||<^J)) 

end 
end 



20 
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EXAMPLE 9 

An experiment was performed to assess the improvement in accuracy of the 
average normal curvature estimation technique over mean curvature estimatioiL 
Synthetic spheres of varying diameter were generated, tessellated using the modified 
marching cubes method, and filtered using two iterations of Algorithm 4.5. Curvature 
estimates were computed based on average normal curvature and mean curvature. The 
reduction in the variance of curvature estimates was used to show the improvement in 
curvature estimation by average normal curvature over mean curvature, as the true 
curvature of a sphere should be constant (zero variance). Figure 54 shows the results of 
this experiment as the percent reduction in the variance of curvature estimates when 
using average normal curvature over mean curvature. There is a minimum of 50% 
percent reduction in variance (increase in reproducibility) when using the average normal 
curvature estimates. This advantage increases with decreased sphere size, as the 
tessellation becomes more coarse. 

Curvature estimates of each type may be made for each triangle in the surface by 
averaging the curvature values for each of the constituent vertices, 

C^.(Q^alC^L^;T,={K,V,,VJ (6.57) 

This provides an elegant method for visualization of the 3D curvature map. 
Each triangle in the surface may be assigned a color (or gray-level) proportional to the 
curvature estimate determined for that triangle. The complete curvature estimation 
technique is described in Algorithm 6.6. 
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Algorithm 6.6 (Three-Dimensional Surface Curvature Estimation) 



for =V^ :K„ (for each vertex) 
for (y,,V,)GadjXK) 

^3 Tj={V,,V^,V^}(V,isa member of triangle r,. ) 

end 

for Tj =To :T^ (for each triangle T. ) 

xj/j =(za X ib )/za x ib (calculate triangle sur&ce normal) 

end 

<l>j =(^J^^y^y)/(m=l)(calculate vertex surface normal) 
end 

for V. = Vq : V„ (for each vertex) 
for V^^adJiV,) 

0,=cos-^((0,-(^J/(l(|.,|l(#.J)| 

= ^a^^^ea (calculate curvature at vertex) 

end 

end 

for Ti =To :T„ (for each triangle) 

Ti={Ya.V,,VJ 

—yCiTa "^^Kft //3 (calculate curvature on triangle) 

end 

Figure 55 shows the result of 3D surface curvature estimation on a malignant 
small pulmonary nodule. The brighter sxufaces correspond to regions of higher 
convexity, whereas the darker siirfaces are regions of higher concavity. The intermediate 
gray-level between white and black is assigned to flat regions. 

Histogram Analysis of Curvature Estimates - With methods to 
estimate the surface curvature of the nodule, we may now explore how such 
data can be analyzed and made suitable as a nodule metric. In particular, 
we would like to devise a curvature metric appropriate for the 
differentiation of benign from malignant nodules. 

It has been documented in the literature that malignant nodules are more bumpy 
(spiculated, lobular) than benign ones [26, 82, 73]. One appropriate goal for a 
curvature-based metric, therefore, would be to characterize the degree of surface 
irregularity present in pulmonary nodules. This may be accomplished using the surface 
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curvature estimation method described above, combined with an analysis of the 
frequency distribution of values measured in each pulmonary nodule. 

Curvature Distribution Model - Consider a perfect sphere of radius r. At each 
point on its surface, the mean curvature has a constant value of 1/r. The curvature of a 
sphere, then, is constant and inversely proportional to its size. 

In our curvature model defined in 3D tessellated space, the average normal 
curvature estimate of a "perfect sphere" asymptotically approaches a single value as the 
density of triangles in the tessellation increases For practical spheres in our tessellated 
model, the average normal curvature will have a range of values whose variance 
decreases with sphere size. Furthemiore, the mean of this distribution will approximate 
the mean curvature of the sphere, or where r is measured in voxels. Figure 56 
shows tiie mean average normal curvature estimate for synthetically-generated spheres 
of up to 10 cm in diameter, sampled on a perfectly isotropic 0.25 mm grid, tessellated 
using Algorithm 4.4, and filtered using two iterations of Algorithm 4.5. The theoretical 
curvature values by sphere diameter are shown using the deished line in the same graph. 
Given the 0.25 mm resolution of the grid, the radius of each sphere in voxels is Adl2. 
Consider the 1 mm sphere. Its radius is 2 voxels, and the corresponding curvature 
estimate is 1/2 = 0.5 radians. Similarly, for the 5 mm sphere, the mdius is 10 voxels 
and the curvature estimate is 0. 1 radians. 

We may begin with a nodule model that is spherical and determine how 
perturbations in the nodule surface (such as those arising from spiculations) affect the 
distribution of the curvature estimates made on the nodule representation. Figure 57 is a 
2D illustration of a nodule with a surface perturbation. In this simple illustration, the 
surface is perturbed to contain a smaller region of constant curvature (circular/spherical). 
It can be seen that in this case, the frequency distribution of curvature estimates will have 
three basic components, a large peak corresponding to the general curvature of the nodule 
and two additional peaks, one containing the curvature of the protrusion (B), and one 
corresponding to the local concavity, C, at the interface between convex curvatures A and 
B. 
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Consider the following illustrative example. The surface of a spherical nodule 
model of diameter dn is perturbed using a local spherical function at a number of points, 
n. The diameter of each perturbation is dp. Given our model, we wovdd expect three 
curvature peaks in the firequency distribution of curvature estimates for the surface, one 
5 each corresponding to the curvature of the original nodule surface, to the curvature of the 
perturbations, and to the concavity introduced at the interface between the large sphere 
and each perturbation. We will consider an example where dny = 40, dp = 8, and n — 26. 

Figure 58 shows the frequency distribution (histogram) of average normal 
curvature estimates in a synthetic sphere which is 40 voxels in diameter. The theoretical 
10 curvature of this sphere is therefore. 

1/r = 2/d = 0.05 radians 

and the observed frequency distribution (mean = 0.046, SD = 0.048) agrees, as expected. 
Next, we may examine the distribution of measured curvature estimates for the spheres 
(^7=8) with which the surface will be perturbed. This histogram (mean = 0.241, SD = 

1 5 0.032) is shown in Figure 59. 

The result of performing surface curvature analysis on the entire perturbed sphere 
model; however, contains more than simply the sxmi of the two characteristic curvatures 
of the large (nodule) and small (perturbations) spheres. Figure 60 shows the frequency 
distribution of curvature estimates for the entire model. In addition to the peaks 

20 corresponding to the nodule surface (0.05 radians) and to the curvature of each 

perturbation (0.25 radians), there is also a concave component at approximately -0.2 
radians. This is, in fact, the manifestation of the expected concavities at the interface 
between the spherical sxirface and each smaller spherical perturbation. 

The location of the two peaks added to the histogram as a resxjlt of surface 

25 perturbation increase in absolute magnitude with decreasing perturbation diameter. 
The convex component increases due to its decreasing radius of curvature, and, 
correspondingly, the magnitude of the (negative) concave component increases as the 
angular difference between the surface and perturbation increases at the point 
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of junctioiL As an example. Figure 61 illustrates another perturbed sphere model 
where the diameter of the perturbations is half that used in the example shown in 
Figure 60. Note that the curvature peaks for the perturbation surface and the nodule- 
perturbation interface have each moved away from Ihe origin. 

Using such a model, several observations can be made. First, it can be seen that 
as the number of surface perturbations increases, the number of added curvature peaks 
increases. Second, the curvature of each of the perturbations is higher than that of the 
underlying nodule surface. Third, the curvature estimates resulting from the interfaces 
of each outward perturbation contribute negative (convex) curvature values to the 
distribution. Each of these observations lead to the conclusion that wdth increased 
numbers of surface perturbations and decreasing diameter of each perturbation, the 
width (and correspondingly, the variance) of the curvature distribution will increase. 
Thus, a measure of the variance of the frequency distribution of curvature estimates 
will provide a quantitative means of characterizing the irregularity of the nodule 
surface. Furthermore, as nodule surface irregularity has been shown to correlate with 
malignancy, this cvurvature variance metric may prove useful in the discrimination of 
benign and malignant nodules. 

One more note must be made on the evaluation of the distribution of curvature 
estimates. When analyzmg the 3D curvature of a nodule which has been segmented to 
remove other structures (e.g. vessels), there may be artifacts, such as artificially convex 
regions, at these points of segmentation. For this reason, when 3D curvature analysis is 
performed, the cxirvature estimates for those regions should be discarded. The 
procedure may be done as follows: (a) segment the pulmonary nodule, removing 
extraneous stmctures; (b)generate the 3D tessellated surface model; (c) compute 3D 
curvature estimates for each vertex and triangle in the model; (d) mask the surface 
model of the nodule vising the dilation of each of the deleted components; and (e) 
compute the frequency distribution of the remaining curvature values. 
This was the method used in the 3D curvatme analysis of all of the nodules in this 
work. The results are described in below. 
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Sensitivity of Curvature Estimates -Curvature estimates and the analysis of 
their distribution are sensitive to tiiree general parameters: (i) image resolution; (ii) 
segmentation method; and (iii) degree of surface smoothing 

The original CT image resolution and the isotropically-resampled image 
5 resolution affect the degree of accuracy with which curvature estimates can be made. 
Improved (finer) resolution allows curvature estimates to be made more locally^ at the 
cost of increased noise. Parameters of the segmentation method, including any . 
intensity thresholds used, affect the curvature estimates by redefining the boimdary of 
the nodule being studied. Finally, the degiee of surface smoothing (in this case, the 
10 nxamber of iterations of Algorithm 4.5 affects the curvature estimates by removing 

high-firequency noise due to spatial quantization. Over-application of this smoothing 
process, however, may result in the elimination of true curvature chamcteristics. 

EXAMPLE 10 

E^eriments were performed to assess the affect of Algorithm 4.5 on the 
1 5 distribution of surface curvature estimates. First, a synthetic sphere of 1 0 nam in 
diameter (40 voxels) was tessellated using Algorithm 4.4. A varying number of 
smoothing iterations were then applied and the variance of the surface curvature 
distribution calculated. Figure 62 shows the results of this experiment. It can be seen 
that in each iteration, the variance of the curvature distribution decreases as each 
20 curvature estimate approaches the trae mean value for the sphere. This observed 

reduction in variance is approxirnately proportional to 7/«2, where n is the number of 
smoothing iterations. 

Malignancy Models -Malignancy models are mathematical tools used to predict 
the likelihood of nodule malignancy based on a variety of parameters. The majority of 
25 these parameters may be calculated feature metrics, such as those described in the 

preceding sections. Malignancy models may also include, however, risk-related patient 
demographic information, such as age, smoking history, and exposure to other known 
limg carcinogens. The actual selection of parameters may, therefore, be based on a 
number of methods, mostly statistical, that maximize the ability of the model to 
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characterize the likelihood of nodule malignancy across a wide range of pulmonary 
nodules. 

Statistical Parameter Selection -Parameters in a model of nodule malignancy 
are generally selected in two steps. The jBrst step that is employed is to include 
5 parameters that would be expected to be different in benign and malignant nodules, 

based on the underlying biology. The second step is to select parameters through the use 
of statistical tests tiiat quantify the degree to which each parameter is associated with 
malignancy. 

Student^s t-test was used to compare the distributions of values for each shape 
10 parameter in malignant and non-malignant nodules. This served two pvirposes. The 
first was to determine which parameters, on there own, were significantly different in 
the two cleisses of nodules. This information is quite usefiil in that it can lead to new 
research into refinement of metrics that appear promising as predictors of malignancy. 
The second goal was to determine those parameters which might be successfiUly 
1 5 combined into a model that would be predictive of malignancy. Once promising 

individual parameters were determined, combinations thereof can be used to generate 
logistic regression models for the classification of nodules. 

The logistic regression model is a type of nonlinear multiple regression which is 
used to model data with categorical (non-continuous outputs). For example, a standard 
20 multiple linear regression model, y^f{xuX2^..,^Xr^ uses a linear combination of the 
input parameters, jc„ to predict the output;/: 

jj; = a + b{K^ + 62^2 + — ^if^w (6.58) 
The regression method, generally based on linear least-squares methods, uses a 
known set, of observation vectors, x, and the corresponding vector of outputs, y, to 
25 optimize the coefiScients of the model, a and 61 . . . An- The range of output values in 
linear regression model is continuous, as the output, y could be any real number. 

In the case where the outputs are categorical, however, logistic regression 
models may be used. These allow the outpvrt to be restricted to a particular range, 
typically between 0 and 1 . For cases where the process being modeled has multiple 

30 discrete values, miiltinomial (or polytomous) logistic regression may be used. 
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However, most frequently, as in the case of nodule malignancy models, the dependent 
variable is binary- valued, either malignant or non-malignant (benign). Therefore, a 
binary logistic regression model is appropriate. 

Logistic regression is based on the logit trajisformation, which is the natural log 
5 of the odds ratio. The odds ratio of an event is a function of the probability, P of that 
event: 

odds = (6.59) 

The logit transformation, then is 

r p 

10 Logistic regression models are based on the assumption that the dependent 

variable (e.g. malignant or benign) is linearly dependent not on the observed 
parameters, xu but on the logarithm of the odds ratio. Therefore, 

( P \ 

Xnipdds) = \n =a 4- b^x^ + &2^2 + —^»^» (6.61) 

\l-P) 

It can be seen, then, that the probability of the event in terms of the observed parameters 
15 is 

P = eKp(a + b^x^ + b^x^ + >-&n^J g2) 
1 +exp(a + b^x^ + 62^2 + ) 



logit {P) = Inipdds)^ 



(6.60) 



which is equivalent to 



P= J—. ^- v: (6.63) 

1 + exp(- [a + b^x^ -f A^jCj + )) 



Thus, if our logistic model is 



20 = 7-7 v: i^'^^) 

1 + exp (- [a + + 62^2 + -K^n )) 

the range of outcomes is between 0 and 1 for any values of the observed parameters, 
and furthermore, the output of the model varies sigmoidally between the two extremes. 

The logistic regression model was used in this work to develop malignancy 
models based on two- and three-dimensional shape and density features. Selection of 
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the individual parameters and performance of the models will be discussed in Section 
below. 

EXAMPLE 11 

In order to quantify the ability of each of the shape metrics studied in 
5 characterizing the dijSerences between benign and malignant small pulmonary nodules, 
a series of tests were performed using a group of in vivo puhnonary nodules, scanned 
using high-resolution CT. The metrics were first compared individually, to assess 
differences in the distributions of each parameter in the two nodule classes. Next, 
logistic regression models were developed based on 2D and 3D metrics, as well as 

10 metrics derived from 3D surface curvature analysis. 

Individual 2D and 3D Metrics - As a first step in assessing the suitability of 
each shape metric for the prediction of malignancy. Student's t-test was applied to test 
the difference in means of the distribution of each parameter for benign and malignant 
nodules. Table 6.2 shows the significance of the difference in means between values of 

1 5 size-variant 3D shape metrics for 22 benign and malignant pxilmonary nodules In 

addition, the means and standard deviations of the distribution of each metric for benign 
and maUgnant nodules is given. As size has been well-documented as associated with 
malignancy, it is not surprising that each one of these metrics significant (p < 0.05), and 
that most are quite significant (p < 0.01). 

20 Table 6.2: Analysis of distribution of size-variant 3D metrics for benign and maUgnant 
nodules 



Parameter 


P Value 


Benign 


Malignant 


Mean 


SD 


Mean 


SD 


height 


0.0015 


3.53413 


0.81352 


6.02052 


2.39888 


surface area 


0.0031 


96.309 


62.419 


339.126 


271.547 


x-extent 


0.0037 


5.7000 


1.88557 


10.6429 


5.26952 


z-extent 


0.0038 


5.05000 


1.39578 


8.64286 


3.81842 


width 


0.0054 


5.16328 


1.63050 


8.21210 


3.00452 
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length 


0.0081 


6.6068 


2.47830 


10.7895 


4.22890 


volume 


0.0087 


75.546 


66.686 


383.776 


410.535 


mass 


0.0096 


11931.1 


12136.2 


57878.3 


61302.5 


y-extent 


0.0178 


6.21667 


2.43645 


9.78571 


4.06568 



Table 6.3: Analysis of distribution of size-variant 2D metrics for benign and malignant 
nodules 



Peirameter 


P Value 


Benign 


Malignant 


Mean 


SD 


Mean 


SD 


lei^;th 


0.0010 


5.73953 


1.94499 


9.96925 


3.20678 


x-extent 


0.0071 


5.46667 


1.83193 


9.71429 


4.90626 


y-extent 


0.0194 


5.43333 


2.02543 


8.42857 


3.53764 


mass 


0.0240 


1086.72 


627.56 


2357.92 


1842.10 


width 


0.0247 


4.57575 


1.18593 


6.70308 


2.98820 


area 


0.0444 


36.1917 


37.4608 


81.4286 


61.6513 


perimeter 


0.2400 


33.7924 


43.8789 


66.6005 


84.7562 



For comparison. Table 6.3 shows the results for the same 22 nodules using 2D 



5 size- variant features computed on the image of maximum cross-sectional area for each 
nodule. Note that the distributions of each of the nine 3D metrics exhibit more 
significant dififerences between benign and malignant nodules than all but two of the five 
2D metrics. 

As the focus of this work is on small pulmonary nodules, such as those that may 
10 be detected in a limg cancer screening program, we may not always expect size to be a 
useful parameter. This is due to the fact that nodules may be detected at any point in 
their growth progression. Therefore, malignancies may frequently be caught while they 
are still small (especially as incidence cases, see Section 5). Furthermore, in those cases 
that are of larger size (i.e. diameter > 3 cm), malignancy is less often in doubt, and 
1 5 biopsy is more straightforward. For these reasons, this study is primarily restricted to 
nodules with diameters smaller than 1 cm. With that in mind, size-invariant metrics 
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should offer a better hope of differentiating small benign and malignant nodules that are 
all of a similar size. Table 6.4 shows the results of t-test analysis of size-invariant 3D 
shape metrics applied to 22 small pulmonary nodules. The metrics are listed in 
ascending order of p value (descending order with respect to significance). 
5 Results for a comparative group of size-invariant 2D metrics are shown in Table 

6.5. When possible, the analogous 2D metrics to those studied in 3D were included 
(e.g. compactness). Here, the difference between 3D and 2D shape metrics is less 
clear. As with the 3D size-invariant metrics, there are relatively few that exhibit high 
significance when used as the sole parameter in differentiating benign from malignant 

1 0 nodules. Therefore, we may turn to another, more appropriate measure of the relative 
benefit of 2D versus 3D nodule shape metrics. Groups of each may be included in a 
logistic regression model and the relative qualities of the models compared. 

Logistic Regression Models - To test the ability of multiple feature metrics to 
characterize differences between benign and malignant nodules, logistic regression 

1 5 models were developed based on metrics whose individual t-tests showed promising 

differences in distribution between the two nodule classes. In addition^ the models were 
restricted to contain parameters that captured uncorrelated information. Models were 
developed using 2D and 3D metrics to compare the abiUty of each of these groups in the 
characterization of nodxile malignancy. Each was based on three features, 

20 The 2D model included LJVR, an aspect ratio based on the ellipse of inertia of 

the nodule cross-section; APR, an alternate measure of 2D nodule compactness; and 
voxel skewness, a measure of the nomiality of the distribution of density within the 
nodule. 

Table 6.4: Analysis of distribution of size-invariant 3D metrics for benign and 
25 malignant nodules 



Parameter 


P Value 


Benign 


Malignant 


Mean 


SD 


Mean 


SID 


normalized volume/surface ratio 


0.0053 


0.42401 


0.01615 


0.38916 


0.03701 


3D compactness 


0.0064 


0.81399 


0.09248 


0.64124 


0.17669 
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volume/surface ratio 


0.0124 


0.69891 


0.14654 


0.96051 


0.30697 


sphericity 


0.1861 


0.47675 


0.16617 


0.37546 


0.15039 


A 

Ad 


0.1888 


6.0125 


14.043 


50.9941 


130.101 




V. 1 J' Of 










Wu 


0.2234 


4.0257 


1 1 .4424 


25.9638 


67,4143 


yaw 






A1 177'? 




A/=; '^ROl 


• 

variance 


U.jO/O 






OUj.U^O 


»>00.\/^o 


roU 


0.3846 


-3.1443 


75.422 


-38.859 


111.360 




0.5462 


0.03285 


0.02325 


0.04076 


0.03711 


WHR 


0.6037 


1.48294 


0,45815 


1.38355 


0.27398 


pitch 


0.6100 


24.6851 


18.9031 


31.1248 


40.2918 


A 


0.6121 


0.00906 


0.00816 


0.00716 


0.00782 


LWR 


0.6519 


1.27539 


0.18937 


1.31137 


0.12034 


YZR 


0.7101 


1 .20822 


0.27464 


1.14632 


0.50272 


standard deviation 


0.7378 


20.2685 


6.0980 


21.6028 


12.6127 


mean 


0.7424 


144.513 


40.7856 


137.408 


57.8951 


LHR 


0.7463 


1.86727 


0.52547 


1.79719 


0.28645 


skewness 


0.7777 


0.07568 


0.32579 


0.11836 


0.32586 


kurtosis 


0.9956 


2.20742 


0.36081 


2.20840 


0.44150 



Table 6.5: Analysis of distribution of size-invariant 2D metrics for benign and 
malignant nodules 



Parameter 


P Value 


Benign 


Malignant 


Mean 


SD 


Mean 


SD 


LWR 


0.0046 


1.24344 


0.17693 


1.58450 


0.32985 


circularity 


0.0313 


0.64667 


0.26691 


0.37888 


0.21562 


skewness 

mm 


0.0852 


-0.626 


0.52973 


-0.1523 


0.65853 
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ZD compaCuicss 


U.Uo /D 




n 074.9 


0 55109 


0 28999 




0.0953 


1.18576 


0.31016 


1.56414 


0.72021 


napr 


0.1722 


0.24347 


0.06839 


0.19826 


0.07285 




0.1994 


0.04062 


0,02536 


0.06981 


0.07872 


XYR 


0.2695 


1.04246 


0.222191 


1.9934 


0.43396 


/A 


0.4483 


4.23139 


11.3991 


0.84959 


0.9959 


w 


0.4o7o 


AO. TX OiC 


1 

iw.D /i 






yaw 


0.4889 

* 


-42.213 


45.7468 


-24.268 


73.6191 


mean 


0.6055 


1 52.847 


43.1481 


141.136 


59.8249 


variance 


0.7178 


481.253 


305.517 


544.265 


502.266 


kurtosis 


0.7660 


0.33915 


25.7179 


3.31157 


0.8487 


standard 
deviation 


0.9083 


20.9205 


6.8338 


20.4526 


12.1223 


A 


0.9466 


0.00166 


0.00134 


0.00160 


0.00298 



The 3D model included 3D compactness, a gross measure of surface 

A 

irregularity; normalized density skew (<^> ^d, a measure of the irregularity of the density 
distribution within the nodule volume; and the variance of tiie average normal surface 
5 curvature metric, a fine measure of surface irregxilarity. 

The pertinent shape metrics were computed for each of 22 small pulmonary 
nodules (7 malignant, 15 non-malignant). Statistical analysis was performed to assess 
the ability of each model to differentiate malignant from non-malignant nodules. The 
results of this experiment are summarized in Table 6.6. 
1 0 Table 6.6; Performance of 2D and 3D metric-based models in characterizing 

malignant and non-malignant small pulmonary nodules 



Nodxile Status 


2D Metrics 


3D Metrics 


Malignant 


Benign 


Malignant 


Benign 
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Malignant (7) 


6 


IFN 


7 


0 


Non-Malignant (1 5) 


1 FP 


14 


IFP 


14 




0.690 


0.778 


P 


<0.0003 


<0.0001 



The model based on 2D metrics was able to characterize much of the difference 
between benign and malignant nodules. The R2 value (used to assess the proportion of 
the data explained by the model) was 0.6900, and the p value (probability that the effect 
5 would be produced by chance) was 0.0003. The model was able to classify all but one 
malignancy (1 false negative (FN))i and one benign nodule (1 false positive (FP)) 
correctly. Therefore, the sensitivity was 86% and the specificity was 93%. The 
positive predictive value was 86%. 

The model beised on 3D metrics fit the experimental data even more closely. 

10 The value of R2 was 0.7776, with p < 0.0001 . In addition, the model was able to 
classify each of the cases correctly, with the exception of one benign nodule (1 FP). 
The sensitivity was 100%, the specificity was 93%, and the positive predictive value 
was 88%. This is an especially good result considering that the penalty in producing a 
false-positive result (unnecessary biopsy) is much less severe than a false-negative one 

1 5 (missed cancer). 

An «-way (leave-one-out) cross-validation approach was also used to evaltiate 
the performance of the model based on 3D metrics, as it may produce a somewhat less- 
biased evaluation of the model- In this experiment, the model was able to correctly 
classify 86% of the nodules correctly, producing 2 FP and 1 FN. The sensitivity was 

20 86%, the specificity was 87%, and the positive predictive value was 75%. To provide a 
more general evaluation, Receiver-Operating Characteristic (ROC) curves were 
generated based on the ability of the model to characterize the entire dataset, as well as 
each of the individual models developed in the cross-validation scheme [27, 2]. The 
area under the ROC curve, Az, for the overall model was 0.9905. The ROC curve for 

25 the complete dataset is shown in Figure 63. The mean area imder the ROC curves 

using the cross-validation scheme was 0.9870, with 95% confidence intervals between 
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0.9822 and 0.9918, while the area under the aggregate curve was 0.9203, TTiis 
aggregate ROC curve is shown in Figure 64. Overall, the model showed considerable 
ability to characterize the cases m the study. 

EXAMPLE 12 

5 Analysis of 3D Surface Curvature - To evaluate the usefulness of the 3D 

sur&ce curvature analysis itself, an experiment was performed using the set of 22 in-vivo 
pulmonary nodules used to evaluate the 2D and 3D feature metrics. The measiire of the 
three-dimensional surface curvature was computed using the techniques described above, 
including the analysis of surface curvature distribution. Each of the nodules studied had 

10 complete 3D information derived from high-resolution (Irom slice thickness) helical CT. 
Curvature values computed from portions of the surface in contact with other structures 
(e.g. vessels) were removed from consideration. Frequency distribution of this measure 
as well as the mean, variance, coefficient of variation, and skewness were detemiined, as 
described hereinabove. The distribution of each of these metrics in the 22 iti-vivo 

15 puhnonary nodules were examined. The results are summarized in Table 6.7. 

Table 6.7: Analysis of distribution of 3D curvature-based metrics for benign and 
malignant nodules 



Parameter 


P Value 


Benign 


Malignant 


Mean 


JSD 


Mean 


SD 


variance 


<0.0001 


0.01084 


0.00395 


0.02274 


0.00499 


standard deviation 


<0.0001 


0.10234 


0.01983 


0.15001 


0.01653 


mean 


0.0003 


0.10349 


0.02682 


0.05095 


0.02419 


range 


0.0006 


0.85084 


0.26098 


1.94865 


0,99224 


skewness 


0.0059 


0.35311 


0.55021 


-0.5373 


0.78756 


kurtosis 


0.0492 


4.40753 


2.40237 


6.99628 


3.29447 


maximum 


0.0999 


0.58676 


0.19492 


0.85648 


0.54787 
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The mean curvature estimate in malignant nodules was 0.0509 ±0.010 (mean db 
sem) vs. 0.1035 ± 0.007 in benign nodules (p <0.0003). The variance of curvature in 
malignant lesions was 0.0227 ± 0.002 vs. 0.0108 ± 0.001 in benign ones (p .0001). The 
skewness of the curvature distribution was -0.537 ± 0.298 vs. 0.353 ± 0.142 in malignant 
5 and benign nodules, respectively (p <0.01). Finally, as a mean-normalized measure of 
curvature distribution, the coefficient of variation (CV) of curvature estimates for 
malignant nodules was 3.823 ± 0.919 vs. 1.095 ± 0.132 in benign nodules (p = 0.0004). 
Thus, each metric of the measure of 3D curvature was useful in differentiating between 
small malignant and benign nodules. 

10 Region Analysis Without Segmentation - Image segmentation is one of the 

most challenging steps in the analysis of small pulmonary nodules. It may be difficult to 
remove connected structures (e.g. vessels) without removing information about the 
surface of the nodule itself. In addition, when studying the change la the nodule over 
time (as in volumetric doubling time estimation), it is of great importance that the 

* 

1 5 segmentations of each image are consistent 

It may be possible to study the change in nodule size and shape over time, without 
performing explicit segmentation, however. This chapter describes techniques for 
comparing three-dimensional nodule regions of interest (ROIs) as two density 
distributions, without segmentation. 

20 Three-Dimensional Region Registration - Prior to comparison of two three- 

dimensional ROIs, they must be well-defined so that each corresponds to the same limg 
region in each study. This is a problem of three-dimensional region registration. Two 
techniques for region registration will be discussed: 1) center of mass (COM) alignment 
and 2) correlation-based region matching. 

25 Center of Mass Alignment -The two 3D ROIs may be aligned using a simple 

procedure based on the computation of the center of mass. Beginning from an initial 
starting point in the 3D image, an iterative search method is used to find the location of 
the COM in each of a number of translated ROIs. In each iteration, the current location of 
the COM, Cc is compared with the actual COM measxired for the current ROI, Cm. The 
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next location of the COM is then replaced with C,n and the ROI recomputed. The 
stopping condition is when the location of the COM does not move by more than some 

small distance € between iterations. Algorithm 7.1 describes the center of mass 
determination technique, where D(Va, is the Euclidean distance function between 
5 voxels Va and P% 

Algorithm 7.1 (Iterative Center of Mass Determination) 

Select an initial location for the COM: C/ = v (jc i, yu ^i) 
Select the desired size of the ROI: S — x,y,z 

AC = oo 
10 Cc=C/ 

whae(AC>e) 

Compute new ROI bounds based on Cc and S 

Cm (X) = /Wooo/'WOOO 

Cm (y) = moio/wooo 
IS Cm (z) = wooi/mooo 

^c^D{Cn,lCc) 

Cc ^ Cm 

end 

The inputs to this procedure, in addition to the 3D image, are an initial starting 
20 point for the search, and the dimensions of the ROI desired. A two-dimensional example 
of the operation of Algorithm 7. 1 is shown in Figure 65. Each image represents one 
iteration of the algorithm. In each frame, the white rectangle illustrates the current ROI, 
in which the white cross shows the computed center of mass, C„|. 
This value then becomes the center of the ROI in the next iteration. The algorithm stops 

25 after 10 iterations, as the COM did not change by more than €= 2 voxels between the 

final two iterations. 

Algorithm 7.1 can be used for the detemiination of two matching ROIs in 
sequential scans of a pulmonary nodule, allowing for comparative analysis. The method 
is quite sensitive to the exact distribution of density, including that corresponding to 
30 structures in the periphery (such as vessels), however. In particular, when a nodule 
exhibits non-uniform growth, the resultant ROIs may not correspond perfectly. 
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Correlation-Based Region Matching - In the method described in the previous 
section, we aligned two images by detemiining a robust center of mass for each and 
aligned these, using the assumption that the COM would not shift significantly between 
scans. An alternative approach to the detennination and alignment of two comparable 
5 ROIs may be defined in terms of image correlation. The two ROIs may be aligned using 
3D correlation as a measure of image correspondence. The method requires selection of 
the two ROIs and specification of the limits on the extent of the search. The second ROI 
is then translated to all locations in the search region and the 3D correlation computed 
between the two ROIs. The translation yielding the greatest value of this match metric is 
10 used to specify the alignment of the ROIs for subsequent analysis. The technique is 

described in Algorithm 7.2, where the T (ROI, i, J, k) specifies the translation of an ROI 
coordinate system byx = hy='j\z — k and the star (jk) operator is the three-dimensional 
correlation described in Equation 7.1 . 

c(/,7, ^)=2ZZi /(^>J^>2)'g(x-/,:>/-7, z-^) (7.1) 

X y z 

15 In each iteration, the correlation of the two ROIs at the current translation of ROIb 

is computed. If the value of this correlation, c, is greater than the previous maximxmi 
value, M, the value and the location of this better match, Z, are recorded. The final result 
of the algorithm is that translation, i, that produces the best alignment between the ROIs. 



20 Algorithm 7.2 (Correlation-Based Region Matching) 

Select two ROIs: ROIa, ROIb 

Select the extent of the search region in each dimension: iS* = {jc, y, z) 

M=0 

for k — ^SziSz {for each translation in search region} 
25 for j = 'Sy : Sy 

for i^Sx iSx 

c(/,y; k) ^ ROIa ROIb 
ifcaj\k)>M 
M^c{iJ, k) 
30 L^ihlk) 

end 

end 

end 
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end 



An example of the result of correlation-based registration is shown for two scans 

of a small pulmonary nodule in Figure 60. When using correlation-based matching, 

changes in the absolute center of mass of the nodule are less likely to affect the 
5 registration, as the overall match is based on direct correspondence between areas of 

intensity in each ROI. For this reason, the correlation-based matching technique was 

chosen for the analysis methods described later in this chapter. 

Three-Dimensional Region Weighting - Since we are analyzing the nodule ROI 

without performing explicit segmentation, vessels and other confounding structures will 
1 0 influence the metrics computed on the region as a whole. In an effort to reduce the 

influence of peripheral structures, we may prefQter the region using a weighting function 

that gives the highest weight to central structures (near the center of mass), and less 

weight to the periphery. In this section we will exanune several weighting functions and 

their effects on our model of the nodule region. 
15 Weighting Functions - Two region weighting functions were considered when 

developing the region analysis technique. These included: (a) rectangular radial window; 

(b)triaugular radial window; and (c) Guassian window. 

In the following discussion, we will assume the ROI to be translated such that the 

center of mass is located at the origin of the coordinate system, which simplifies the 
20 mathematical notation. With that in mind, the rectangular radial window can be defined 

as follows: 



In the one-dimensional case, this degenerates to the familiar recto 0 function, 
where the width of the rectangle is 2r. In two dimensions, the window is circular with 
25 radius r. Similarly, in the three-dimensional case the window is spherical with radius r. 
This function, while easy to implement has several disadvantages. First, it imposes a hard 
cut-off point, giving zero weight to any voxels outside of r. More important, there is no 




(7.2) 
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scaling given to any of fhe interior region. This, in effect, imposes a hard segmentation 
based on distance &om the center of mass. 

An improvement over the rectangular window would be one that provides a 
gradual attenxiation of the voxel intensities as they appear further from the COM. Such a 
5 window would allow the weighting function to retain the majority of the density 

information near the center of ROI and much less near the periphery. This lead to the 
choice of a Gatissian weighting function. 

The one-dimensional Gaussian may be specified ia terms of it's standard 
deviation, a. The intensity at a point relative to the mean, /i, is defined as 

10 a)g(x; = — ^ e\x-^if /2a^ (7.3) 

In two dimensions, the distribution may have a different standard deviation for 
each dimension (resulting in an elliptical Gaussian), but as we wish to obtain a symmetric 
weighting function, we will consider the circular case where 

15 The corresponding two-dimensional Gaussian function may be expressed as 

y) = e-K^-''/7/2a^ (7.4) 

Zttc 

where the origin is located at ( ii^^fXy). Similarly, we may define a symmetric (spherical) 
three-4imensional Gaussian weighting function as 

rx. y, z; = e -K'-^)^^-''.)'^-'^)^!'-^ (7.5) 

20 Whereas we will set the mean (origin) of the Gaussian weighting fimction to be the center 
of mass of the ROI to be weighted, the choice of an appropriate value for the standard 
deviation, a, is less obvious. The overall goal of the weighting is to provide a 
representation of the original ROI with structures at the periphery smoothly attenuated, so 
that they contribute less to the overall density distribution. Thus, the choice of values for 

25 a must be a function of two parameters: the size of the ROI, and the degree to which the 
periphery should be attenuated. In addition, the weighting function used for an ROI 
containing a nodule in one CT scan shoxild be the same as that used to weight the ROI in 
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a subsequent scan. A discussion of the appropriate mechanism for selection of values for 
G will be given in the context of the doubling time estimation problem in Section 7.3.1. 

Moment-Based Analysis - The method of moments may be employed to study 
the density distribution of the weighted ROI. Moments were discussed in Section 6. In 
5 this context, we are using densitometric moments as descriptors of the weighted region 
containing the nodule. First and foremost is the zeroth-order moment: 

mooo^Yi Z S ^i^^y^^) (7.6) 

jc=0 y^O r=0 

where v(x, y, z) is a continuous, real-valued function describing the CT attenuation at a 

given location in 3 -space. Recall that when v (x, y, z) is binary-valued, as in a segmented 
10 ROI, the zeroth-order moment corresponds to the volume of the nodule region. Li this 

case, however, it is the sum of the attenuation at each voxel in the ROI. More 

specifically, if we assxmie that CT attenuation corresponds to object density at a particular 

location, /wooo niay be thought of as a measure of region mass. 

In the current application, however, we are analyzing a weighted nodule ROI. 
15 The goal, then, is to determine a moment-based descriptor that can be used to 

characterize the nodule or change in the nodule over time. We will, therefore, define the 

mean density of the weighted nodule ROI as 

This is simply the "mass" of the region divided by the volume of the ROI, 
20 analogous to physical density measurements. Using this metric, we are able to study the 
change in size of the nodule over time. 

Doubling Time Estimation -Doubling time estimates based on the change in 
segmented nodxde volume were discussed in Section 5. An alternative approach to 
doubling time estimation not requiring explicit segmentation would be to quantify the 
25 change in mean density of the ROI as a surrogate for change in nodule volimie. 

Consider two scans of an growing object at two times within ROIs of identical 
size. At baseline, the volume of the object is Ki, and in the second scan it is V2. If we 
consider the case v^diere the object is of uniform density and that the background is zero- 
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valued (empty), the ratio of the mean densities of these regions is exactly equal to the 
ratios of their volumes: 

md^ =J^,md^ = J!L=>i!!^ = ^ (7.8) 

This mathematical relationship suggests that it may be possible to use the mean 
5 density metric as a proxy for nodule volume in computing doubling time estimates. Of 
great importance, obviously, is the selection of two corresponding ROIs with identical 
dimensions. There are two additional considerations, however: (a) the nodule may not be 
of uniform density; and (b) non-nodule objects may be included in the ROI. 

The first consideration supposes that if the nodule changes in density but not in 

1 0 "size," the relationship between mean density and segmented volume no longer applies. 
This is true. However, it may be argued that an increase in measured nodule density, for 
eqtdvalent nodule "size," may be a better indicator of increased numbers of cells in the 
region. Furthermore, the removal of the intensity-based thresholding step in traditional 
nodule segmentation allows us to consider the entire nodule volimie in every case, 

1 5 eliminating the complexity and algorithmic (or radiologist) bias in segmentation. 

The second consideration, pertaining to the inclusion of non-nodule regions of 
significant attenuation in the ROI is also valid. However, most of the non-nodule 
structures in the ROI correspond to vasculature or snudl airways. In each case, the size of 
such structures is likely to remain constant between scans, unless the nodule (as in the 

20 case of some malignant tumors) is responsible for angiogenesis in the region. In this 
case, the additional density measured in the ROI would contribute to a more aggressive 
estimate of growth, appropriate in such cases. Furthermore, the contribution of structures 
at the periphery of the ROI will be limited as the region has been pre-weighted using a 
Gaussian weighting function. This will help achieve the goal of keeping the density 

25 distribution of the nodule itself as the chief contributor to the density distribution of the 
entire ROI. 

Methods - A complete method for nodule doubling time estimation without 
explicit nodule segmentation has been developed. The goal is to use the mean density of 
the Gaussian-weighted nodule ROI as a proxy for the true nodule volume. In this way, 
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we can compare the mean densities of two appropriately weighted nodule ROIs taken at 
two different times and determine an estimate of the rate of growth of the nodule. A key 
parameter in this procedure is the standard deviation, a, used in each of the ROI 
weighting functions. It is desirable to have the Gaussian taper toward zero near the 
5 periphery of the ROI so that additional, non-nodule structures (e.g. vessels) present in the 
ROI will have less of a contribution toward the density calculations. This would suggest 
a relatively small value of o might be appropriate. However, we do not wish to over- 
attenuate the true nodule region. Furthermore, as our growth estimation technique is 
based on measures of mean density in the region, we do not wish to choose a value of a 

1 0 small enough such that the mean density metric will become "saturated" upon subsequent 
nodide growth. This might otherwise be possible as we are using the same weighting 
function for the ROI at each time. For example, if the first ROI is weighted such that the 
mean density is quite high (near the effective maximum of the metric), nodtile growth 
that appears in the second ROI may be missed if the radius of the Gaussian is too small. 

15 Furthermore, increases in mean attenuation near the center of the nodule in the second 
ROI may not be accurately measured in this case, as there might be insufficient dynamic 
range in the weighted density scale. 

For these reasons, an iterative technique is employed to determine the value of a 
for the Gaussian weighting function such that the mean density in the first ROI is fixed at 

20 a pre-determined point. In this way, a sufficient amoxmt of dynamic range in the mean 
density metric is preserved for measurement of growth in the second ROI Algorithm 73 
describes the method for selection of the appropriate value of a for a given standardized 

mean density, md^, and value of e, the stopping criterion. 

Algorithm 7.3 (Iterative Sigma Determination) 

25 Select an ROI, ROI, and a standardized mean density value, rndg 

o^b =crc = 0 

J = 0.1 md, =(? 

m<ic= 0 

while done ^1 
30 Generate Gaussian weighting fimction using 

WROI, (x, z) = {x,y, z) • ROI(x,y, z) 
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x=0 y=0 r=0 

if |5|>0 



if |5|<e 

done =1 
else 



10 else 



end 

if|5|<e 

done = 1 
else 

s = s/2 

15 cr^=<7^+^ _ 

end 
end 
end 

At the teimination of the algorithm, the value of a is the appropriate standard 
20 deviation for the Gaussian weighting function, such that it will produce the desired mean 
density, mds. 

■ 

Using the mean density of the weighted nodule ROIs as a surrogate for the nodule 
volume, it is possible to estimate the growth rate of a nodule in sequential CT scans 
without performing explicit segmentation. Given two 3D nodule ROIs, ROIj and ROI2, 

25 acqiiired at // and t2, respectively, we may use the first to compute the appropriate 

weighting function for the pair, and subsequently produce flie Graussian-weighted regions 
of interest, WROIi and WROI2. In each of the Gaussian-weighted 3D regions, we then 
compute the mean voxel densities, mdj and md2> Given these mean density measures and 
the time between scans. At, we may compute a region-based doubling time estimate, 

30 DTr as 

Dr,= ^^^'^ (7.9) 
In(md2 1 md^ ) 
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This is the analogous expression to the calculation of nodule doubling time based 
on segmented volumetric data (Equation 5.24). The entire doubling time estimation 
technique is described in Algorithm 7.4. 



Algorithm 7.4 doubling Time Estimation Witliout Explicit Segmentation) 

Select two ROIs, ROIj, ROIz at ti and tz respectively 
Align (ROIj, ROI^ using Algorithm 7.2 
Determine g based on ROIj using Algorithm 7.3 
Generate Gaussian weighting function w using a 
W ROIi (jc, z) = >v(x, y, z) • ROI^ {x, y, z) 

WROhix^y.z) = w{x,y,z)^ROI^{x,y,z) 

M-1 

md^ =(Z WROI,(x,y,z))IMN^D 

, M-1 N-l Z-l 
md,^(X Y Y WROI^ix,y,z))IM^N'L) 
x=0 y=0 z=0 

1«2-A^ 



It is important to note that this algorithm was designed chiefly to operate on well- 
circumscribed and vascularized nodules. Nodules with sufficient proximity or 
connectivity to the pleural sinf ace (pleural tail, and juxtapleural nodules) present 
significant ROI registration problems- In particular, both the COM-based and 
correlation-based aligmnent techniques are quite sensitive to ROI selection when a 
significant portion of the pleural wall is included. 

EXAMPLE 12 

The techniques for three-dimensional region analysis without explicit nodule 
segmentation were compared with those based on segmentation. The following is a 
summary of the results. 

Doubling Time Estimation - A study was performed using eleven pulmonary 
nodules for which two high-resolution studies were available. Five were malignant and 
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six benign, either biopsy-proven or exhibiting no growth over two years. Volumetric 
doubling times (DT) were computed based on three-dimensional segmentations of each 
nodule at each time. Doubling time estimates (DTi^ were also computed using the 
change in mean density of the registered ROI, using the techniques described above. The 
results are shown in Table 7.1. 

Table 7.1 : Comparison of in vivo doubling time estimates with and without explicit 
segmentation 





J-flR 




OUXLLlo 


1 


193.5 


104.4 


Malignant 


2 


46.70 


51.09 


Malignant 


3 


182.8 


73.42 


Malignant 


4 


341.4 


177.4 


Malignant 


5 


224.2 


87.04 


Malignant 


6 


4557 


825.5 


Benign 


7 


580.2 


2026 


Benign 


8 


-2894 


-1571 


Benign 


9 


1167 


395.6 


2YNC 


10 


7162 


3335 


Benign 


11 


801.8 


845.9 


Benign 



Comparative statistics were performed to assess the degree of relationship 
between the two types of doubling time estimate. As it was difficult to make the 
assumption both DT and DTr would be Gaussian-distributed variables, the standard 
Pearson correlation was substituted with the Spearman rank correlation. A non- 
parametric measure of correlation, it is based on the difference in rank of corresponding 
values in each variable. If we define Ri to be the rank of one variable in its distribution 
and iS{ to be the rank of the second variable in its distribution, tiie Spearman rank-order 
correlation coefficient, is defined as 



(7.10) 
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This expression can be simplified if we define the sum squared difference of ranks, D as 

^ = E Wi-'S,)^ (7.11) 



'•.^I-TtFtT (7.12) 



The expression for is then simply 

6D 

N'-N 

5 The value of was computed to show the relationship of volumetrically determined 
doubling times to those detennined from change in mean density. In this experiment, 

was equal to 0-9091, P> |r,| = O.OOOL This indicates a very high correlation between the 

values of DT and DTr. Thus, the estimation of nodide doubling time without expUcit 
segmentation may be an effective measure in place of the segmentation-based method. 

10 An examination of the DTr values alone as a predictor of nodule malignancy 

was also encouraging. The density-estimated doubling time values were significantly 
different for mahgnant nodules, 197.73 ± 47.1 (mean =1= sem), vs. benign ones, 2860.7 i 
1059.6 (p < 0.01, using Wilcoxon rank siuns). Again, these techniques for the 
estimation of nodule growth without explicit segmentation appear promising in the 

1 5 differentiation of benign from malignant lesions. 

Thus, while there have been described what presently believed to be the 
preferred embodiment of the invention, those skilled in the art will realize that changes 
and modifications may be made thereto without departing from the spirit of the 
invention, and it is intended to claim all such changes and modifications as fall within 

20 the true scope of the invention. 
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We claim: 

1. A method of providing a three-dimensional image representation of an 
object, said object capable of detectably attenuating a signal resulting from scaiming, 

5 comprising: 

a. scaiming said object with a scaiming-medium capable of producing 
a signal corresponding to a three-dimensional image of said object; and 

b. obtaining a signal corresponding to said three-dimensional image 
of said object; and 

10 c. supersampling said three-dimensional image signal to a degree 

greater than the highest spatial frequency of said signal. 

2. A method as in Claim 1 wherein said scanning-medium is selected from 
the group consisting of x-ray, magnetic field, ultrasoimd, laser, electrical current, visible 

15 light, ultraviolet light, infrared light, and radio frequency. 

3. A method as in Claim 1 wherein said scanning-mediiom is x-ray. 

4. A method as in Claim 1 wherein said supersampling comprises: 

20 a. imposing three-dimensional units of image representation on said 

three-dimensional image signal, and 

b. increasing the number of three-dimensional image representation 
units to reduce contribution of error per unit. 

25 5. A method as in Claim 1 wherein said supersampling is conducted to 

produce three-dimensional isotropic space. 

6. A method as in Claim 1 wherein said object is located in a host body. 

30 7. A method as in Claim 6 wherein said object is a tissue mass. 
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8, A method as in Claim 6 wherein said object is a growth. 



10 



9. A method as in Claim 8 wherein said growth is a tmnor. 



10. A method as in Claim 6 wherein said object is a puhnonary nodule. 



11. A method as in Claim 1 0 wherein said nodule is not greater than 1 0 mm 
in diameter. 

12. A method as in Claim 4 which further comprises segmenting said object 
from other structures found proximal to said object and scanned with said object in step 
"la." 

15 13. A method as in Claim 12 wherein said segmenting comprises subtracting 

attached structures and other objects from said object. 

14. A method as in Claim 12 wherein said segmenting comprises subtracting 
said object from adjacent structures and other objects. 

20 

15. A method as in Claim 1 2 wherein said segmenting fiuHier comprises 
restoring volumetric and surface characteristics of said object. 

16. An article of manufacture for providing a three-dimensional image 

25 representation of an object, said object capable of detectably attenuating a signal resulting 

from scaiming, comprising: 

a machine readable-medium containing one or more programs which 

when executed implement the steps of: 

a. scanning said object with a scaiming-medium capable of producing 
30 a signal corresponding to a three-dimensional image of said object; and 
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b. obtaixiing a sig^al corresponding to said three-dimensional image 
of said object; and 

c. supersampling said three-dim^isional image signal to a degree 
greater than the highest spatial frequency of said signal. 

5 

17. An article of manufacture as in Claim 16 wherein said scaiming-medium 
is selected from the group consisting of x-ray, magnetic field, ultrasound, laser, electrical 
current, visible light, ultraviolet light, infrared light, and radio frequency. 

10 18, An article of manufacture as in Claim 1 7 v^herein said scaiming-medium 

is x-ray. 

19. All article of manufacture as in Claim 16 wherein said supersampling 
comprises: 

15 a. imposing three-dimensional imits of image representation on said 

three-dimensional image signal, and 

b. increasing the mmiber of Ihree-dimensional image representation 
units to reduce contribution of error per unit. 

20 20. A method as in Claim 1 6 wherein said supersampling is conducted to 

produce three-dimensional isotropic space. 

21. An article of manufacture as in Claim 16 wherein said object is located in 
a host body. 

25 

22. An article of manufacture as in Claim 21 wherein said object is a tissue 

mass. 



30 



23 



An article of manufacture as in Claim 21 wherein said object is a growtia. 
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24. An article of manufacture as in Claim 23 wherein said growth is a tumor. 

25. An article of manufacture as in Claim 21 wherein said object is a 
pulmonary nodule. 

26. An article of manufacture as in Claim 25 wherein said nodule is not 
greater than 1 0 mm in diameter. 

27. An article of manufacture as in Claim 20 which further comprises 
segmenting said object firom other structures found proximal said object and scanned with 
said object in step "16a". 

28. An article of manufacture as in Claim 27 wherein said segmenting 
comprises subtracting attached structures and other objects from said object. 

29. An article of manufacture as in Claim 27 wherein said segmenting 
comprises subtracting said object from adjacent stmctures and other objects. 

30. An article of manufacture as in Claim 27 wherein said segmenting 
further comprises restoring volumetric and surface characteristics of said object. 

31. A system for providing a three-dimensional image representation of an 
object, said object capable of detectably attenuating a signal resulting from scanning, 
comprising: 

a. a scanner for scanning an object with a scanning-medium capable 
of producing a signal corresponding to a three-dimensional image of an object; 

b. a receiver which retrieves said signal for processing; and 

c. a processor configured to receive said signal from said receiver and 

supersample said three-dimensional image signal to a degree greater than the 

highest spatial frequency of said signal. 
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32. A system as in Claim 31 wherein said scanning-medium is selected JQrom 
the group consisting of x-ray, magnetic field, ultrasound, laser, electrical current, visible 
light, ultraviolet light, infrared light, and radio frequency. 



33. A system as in Claim 32 wherein said scaoning-medivim is x-ray. 



34. A system as in Claim 3 1 wherein said supersampling of said processor 
comprises: 

a. imposing imits of image representation on said three-dimensional 
image signal, and 

b. increasing the number of image representation units to reduce 
contribution of error per unit whereby estimation of object boundary is improved. 



35. A system as in Claim 3 1 wherein said supersampling is conducted to 
produce three-dimensional isotropic space. 



36. A system as in Claim 3 1 wherein said object is located in a host body. 



37. A system as in Claim 36 wherein said object is a tissue mass. 



38. A system as in Claim 36 wherein said object is a growth. 



39. A system as in Claim 38 wherein said growth is a tumor. 



40. A system as in Claim 36 wherein said object is a puhnonary nodule. 



41 . A system as in Claim 40 wherein said nodule is not greater than 1 0 mm in 
diameter. 



42. A system as in Claim 34 wherein said processor is configured to segment 
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said object from other structures found proximal said object and scanned with said object. 

43. A system as in Claim 42 wherein said segmenting subtracting attached 
structures and other objects from said object. 

44. A system as in Claim 42 wherein said segmenting subtracting said 
object from adjacent structures and other objects. 



45. A system as in Claim 42 wherein said segmenting further comprises 
restoring volumetric and surface characteristics of said object. 



46. A method for three-dimensional segmentation of a region of interest in a 
host body to differentiate an object in said region, comprising: 

a. scanning said region of interest with a scanning-medium capable 
of providing a signal corresponding to a three-dimensional representation of said 
region; 

b. imposing three-dimensional units of image representation on said 
three-dimensional signal; 

c. identifying image representation units as having object signal, 
background signal, and as boxmdary units having both object and background 
signal; 

d. thresholding to separate background from object image; and 

e. conducting three-dimensional connected component analysis to 
identify as images those objects contiguoiis in three-dimensional space. 



47. A method as in Claim 46 wherein said scanning-medium is selected from 
the group consisting of x-ray, magnetic field, ultrasound, laser, electrical current, visible 
light, ultraviolet light, infrared light, and radio frequency. 

48. A method as in Claim 46 wherein said scanning-medium is x-ray. 
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49. A method as in Claim 46 wherein said thresholding comprises iterative 
threshold selection. 

50. A method as in Claim 46 wherein said three-dimensional connected 
component analysis comprises connected component labeling in three dimensions. 

51. A method as in Claim 50 wherein said coimected component labeling is 
selected from one of (i) recursive connected component labeling, and (ii) iterative 
connected component labeling. 

52. A method as in Claim 50 wherein said connected component analysis 
further comprises selecting objects based on said connected component labeling. 

53 . A method as in Claim 50 which further comprises supersampling in three 
dimensions said signal corresponding to a three-dimensional representation of said 
region. 

54. A method as in Claim 50 which further comprises three-dimensional 
morphologically filtering images identified as a result of steps a-e. 

55. A method as in Claim 54 wherein said filtering comprises three- 
dimensional opening of said images. 

56. A method as in Claim 55 wherein said opening comprises erosion 
followed by dilation of said images in threedimensions. 

57. A method as in Claim 54 which further comprises closing of a three- 
dimensional region of backgroimd signal units resulting from steps a-e. 
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58. A method as in Claim 57 wherein said closing comprises a dilation 
followed by an erosion of said background region in three dimensions. 

59. A method as in Claim 55 which further comprises three-dimenjsionally 
5 regrowing images to approximate original surface characteristics of said images. 

60. A method as in Claim 59 wherein said regrowing comprises iterative 
constrained dilation in three dimensions. 

10 61. A method as in Claim 60 wherein the entire morphological filtering 

process is as follows: 

begin with an initial binary image / 

*^(lQSd)® Sd{PcrfoTtn opening using a spherical kernel Sd of diameter d] 
Perform connected component analysis, keeping the component of interest 
1 5 while > = 2 {Number of useful dilations} 

J = J ® jS^ {Perform dilation using a spherical kernel of diameter s) 

J— J A I {Perform a voxel-by- voxel logical AND } 

s = s/p {where is a constant greater than 1 .0 (does not have to be 
an integer)} 

20 end 

62. A method as in Claim 50 which further comprises eliminating structures 
adjacent to images of interest 

25 63. A method as in Claim 62 wherein said image of interest is a nodule and 

said structure comprises thoracic components. 

64. A method as in Claim 63 wherein said nodule is a juxtapleural nodule 
and said structure comprises pleural surface and thoracic wall. 

30 

65. A method as in Claim 62 wherein said eliminating comprises: 
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i. determining in three dimensions angles describing 
orientation of the surface of structure to be eliminated; 

ii. based on the angles found in (i) performing an opening 
operation to detect a majority of said structure to the exclusion of 
said image of interest; and 

iii. three-dimensionally subtracting said structure from said 
image of interest 

66. A method as in Claim 65 wherein said determining of step (i) comprises 
conducting three-dimensional moment analysis. 

67. A method as in Claim 65 which further comprises morphologically 
filtering in three dimensions said image of interest. 

68. A method as in Claim 65 wherein said eliminating is described as 
follows: 

begui with an initial binary image /, 

using moments, determine the orientation of the pleural surface, 
generate a disk-shaped kernel £>, oriented parallel to the pleural surface 

J—(I&D)®D {Perform morphological opening using D} 
K — I'J {Perform image subtraction} 
Continue with iterative morphological filtering. 

69. A method as in Claim 50 which further comprises generating a smoothed 
surface representation of an image of interest 

70 . An article of manufacture for three-dimensional segmentation of a region 
of interest in a host body to differentiate an object in said region, comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: 

a. scanning said region of interest with a scanning-medium capable 
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of providing a signal corresponding to a three-dimensional representation of said 
region; 

b. imposing three-dimensional units of image representation on said 
three-dimensional signal; 
5 c. identifying image representation units as having object signal, 

backgroimd signal, and as boundary units having both object and backgromxd 
signal; 

d. thresholding to separate background from object image; and 

e. conducting three-dimensional connected component analysis to 
10 identify as images those objects contiguous in three-dimensional space. 



15 



71 . An article of manufacture as in Claim 70 wherein said scaxming-mediimi 
is selected from the group consisting of x-ray, magnetic field, ultrasound, laser, electrical 
current, visible light, ultraviolet light, infrared light, and radio frequency. 



72. An article of manufacture as in Claim 70 wherein said scanning-medium 
is x-ray. 



73. An article of manufacture as in Claim 70 wherein said thresholding 
20 comprises iterative threshold selection. 



74. An article of manufacture as in Claim 70 wherein said three-dimensional 
connected component analysis comprises connected component labeling in three 
dimensions. 

75. An article of manufacture as in Claim 74 wherein said connected 
component labeling is selected from one of (i) recursive connected component labeling, 
and (ii) iterative connected component labeling. 



30 



76. An article of manufacture as in Claim 74 wherein said coimected 
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component analysis further comprises selecting objects based on said connected 
component labeling. 

77. An article of manufacture as in Claim 70 which further comprises 

5 supersampling in three dimensions said signal corresponding to a three-dimensional 
representation of said region. 

78. An article of manufacture as in Claim 70 which further comprises three- 
dimensional morphologically filtering images identified as a result of steps a-e. 

10 

79. An article of manufacture as in Claim 78 wherein said filtering comprises 
three-dimensional opening of said images. 

80. An article of manufacture as in Claim 79 wherein said opening comprises 
1 S erosion followed by dilation of said images in three dimensions. 

8 1 . An article of manufacture as in Claim 78 which further comprises closing 
of a three-dimensional region of background signal units resulting from steps a-e. 

82. An article of manufacture as in Claim 81 wherein said closing comprises 
a dilation followed by an erosion of said background region in three dimensions. 

83. An article of manufacture as in Claim 79 which further comprises three- 
dimensionally regrowing images to approximate original surface characteristics of said 
images. 

84. An article of manufacture as in Claim 83 wherein said regrowing 
comprises iterative constrained dilation in three dimensions. 

30 85. An article of manufacture as in Claim 84 wherein the entire 
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morphological fatering process is as follows: 

begin with an initial binary image / 

J^(lQSd) © Sd {Perform opening using a spherical kemel Sd of diameter d} 
Perform connected component analysis, keeping the component of interest 
5 while s>-2 {Number of useful dilations} 

J = J © jS^ {Perform dilation using a spherical kemel of diameter s} 

J- J A I {Perform a voxel-by-voxel logical AND} 

s ~ s/p {where "j?" is a constant greater than 1 ,0 (does not have to be 
an integer)} 

10 end 

86. An article of manufacture as in Claim 70 which further comprises 
ehminating structures adjacent to images of interest 

1 S 87. An article of manufacture as in Claim 86 wherein said image of interest 

is a nodule and said structure comprises thoracic components. 

88. An article of manufacture as in Claim 87 wherein said nodule is a 
juxtapleural nodule and said stmcture comprises pleural surface and thoracic wall. 

20 

89. An article of manufacture as in Claim 86 wherein ssiid eliminating 
comprises: 

i. determining in three dimensions angles describing 
orientation of the surface of stmcture to be eliminated; 
25 ii. based on the angles found in (i) performing an opening 

operation to detect a majority of said stmcture to the exclusion of 
said image of interest; and 

iii. three-dimensionally subtracting said stmcture from said 
image of interest. 

30 

90. An article of manufacture as in Claim 89 wherein said determining of 
step (i) comprises conducting three-dimensional moment analysis. 
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91 . An article of manufacture as in Claim 89 which further comprises 
morphologically filtering in three dimensions said image of interest, 

5 92. An article of manufacture as in Claim 89 wherein said eliminating is 

described as follows: 

begin with an initial binary image /, 

using moments, determine the orientation of the pleural surface, 
generate a disk-shaped kernel Z), oriented parallel to liie pleural surface 

10 J=(/©D)©jD {Perform morphological opening using D} 

K-I-J {Perform image subtraction} 
Continue with iterative morphological filtering. 

93. An article of manvifacture as in Claim 70 which finther comprises 
1 5 generating a smoothed surface representation of an image of interest. 

94. A system for three-dimensional segmentation of a region of interest in a 
host body to differentiate an object in said region, comprising: 

a. a scanner which scans said region of interest with a scanning- 
20 medixma capable of providing a signed corresponding to a three-dimensional 

representation of said region; and 

b. a processor configured to 

(i) impose three-dimensional units of image representation on said 
three-dimensional signal; 
25 (ii) identify image representation units as having object signal, 

background signal, and as boundary units having both object and 
background signal; 

(iii) threshold to separate background firom object image; and 

(iv) conduct three-dimensional connected component analysis to 
30 identify as images those objects contiguous in three-dimensional 

space. 
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95. A system as in Claim 94 wherein said scaiming-medium is selected from 
the group consisting of x-ray, magnetic field, ultrasomid, laser, electrical current, visible 
light, ultraviolet light, infirared light, and radio frequency. 

5 96. A system as in Claim 94 wherein said scanning-mediimi is x-ray. 

97. A system as in Claim 94 wherein said threshold comprises iterative 
threshold selection. 

10 98, A system as in Claim 94 wherein said three-dimensional connected 

component analysis comprises connected component labeling in three dimensions. 

99. A system as in Claim 98 wherein said connected component labeling is 
selected from one of (i) recursive connected component labeling, and (ii) iterative 

1 5 connected component labeling. 

100. A system as in Claim 98 wherein said connected component analysis 
frirther comprises selecting objects based on said connected component l£ibeling. 

20 1 0 1 . A system as in Claim 94 wherein said processor is ftirther configured to 

supersampling in three dimensions said signal corresponding to a three-dimensional 
representation of said region. 

102. A system as in Claim 94 wherein said processor is frirther configured to 
25 three-dimensional morphologically filter said images. 

103. A system as in Claim 102 wherein said processor filters said images by 
three-dimensional opening of said images. 

30 104. A system as in Claim 103 wherein said opening comprises erosion 
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followed by dilation of said unages in three dimensions. 

1 05. A system as in Claim 1 02 wherein said processor is further configured to 
close a three-dimensional region of background signal units. 

5 

106. A system as in Claim 105 wherein said processor closes said background 
region by a dilation followed by an erosion of said background region ia three 
dimensions. 

10 1 07. A system as in Claim 1 03 wherein said processor is further configured to 

three-dimensionally regrow images to approximate original surface characteristics of said 
images. 

108. A system as in Claim 107 wherein said processor regrows said images by 
1 5 iterative constrained dilation in three dimensions 

109. A system as in Claim 108 wherein said processor is configured to 
morphological ly filtCT as follows: 

begin with an initial binary image / 

20 J='(jQSd)® Sa {Perform opening using a spherical kemel S^t of diameter d} 

Perform connected component analysis, keeping the component of interest 
while 5 > = 2 {Number of useful dilations} 

J^J®S^ {Perform dilation using a spherical kemel of diameter s} 

J- J A I {Perform a voxel-by-voxel logical AND} 
25 5 = sip {where "p" is a constant greater than 1 .0 (does not have to be 

an integer)} 

end 

110. A system as in Claim 94 wherein said processor is fijrther configured to 
30 eliminate structures adjacent to images of interest. 



111. A system as in Claim 110 wherein said image of interest is a nodule and 
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said structure comprises thoracic components. 

112. A system as in Claim 111 wherein said nodule is a juxtapleural nodxile 
and said structure comprises pleural surface and thoracic wall. 

5 

113. A system as in Claim 110 wherein said processor is configured to 
eliminate by: 

i. determining in three dimensions angles describing 
orientation of the surface of structure to be eliminated; 
10 ii. based on the angles fo\md in (i) performing an opening 

operation to detect a majority of said structure to the exclusion of 
said image of interest; and 

iii. three-dimensionally subtracting said structure from said 
image of interest. 

15 

114. A system as in Claim 113 wherein said determining further comprises 
conducting three-dimensional moment analysis. 

115. A system as in Claim 1 13 wherein said processor is further configured to 
20 morphologically filter in three dimensions said image of interest 

116. A system as in Claim 113 wherein said processor eliminates structures as 
follows: 

begin with an initial binary image /, 
25 using moments, determine the orientation of the pleural surface, 

generate a disk-shaped kemel D, oriented parallel to the pleural surface 

J— (/ 0 Z) ) ® D {Perform morphological opening using D} 
K^I'J {Perform image subtraction} 
Continue with iterative morphological filtering. 

30 



117 



A system as in Claim 94 wherein said processor is fiirther configured to 
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10 



generate a smoothed surface representation of an image of interest. 

118. A method of generating a smoothed surface representation of a three- 
dimensional voxel representation of a three-dunensional image comprising: 

a. segmenting a three-dimensional image to provide a segmented 
voxelated three-dimensional signal; 

b. modified tessellating said segmented 3D image to provide a 3D 
triangulated surface representation; and 

c. filtering said 3D representation resulting fi-om (b) to smooth said 
3D triangulated siirface representation. 



119. A method as in Claim 118 wherein said modified tessellating is as follows: 



for all N xy^z — v(jc...jc + 1, y,„y + 1, z...z + 1) 
compute octant o{x, y, z) 
index polygon look-up table using o(x, y, z) 
15 add triangles corresponding to o(x, y, z\ offset by (x, y, z). 

end 

120. A method as in Claim 118 wherein said filtering comprises signal 
processing as follows: 

for A:= 1 : n 

20 for all Vi V {for each vertex} 

S,(x,j;,z)<-0 

for all e adfiVi) {for each vertex adjacent to Vi} 

S,(x,y,z)<r-S,'hVj 

end 

25 V, <r- (1 ^a)V, +aCS, I \ adj(y,) \) 

end 
end 

121 . An article of manufacture for generating a smoothed surface 

30 representation of a three-dimensional voxel representation of a three-dimensional image 
comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: 

a. segmenting a three-dimensional image to provide a segmented 
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voxelated three-dimensional signal; 

b. modified tessellating said segmented 3D image to provide a 3D 
triangiilated surface representation; and 

c. filtering said 3D representation resulting fi:om (b) to smooth said 
3D triangulated surface representation. 



122. An article of manufacture as in Claim 121 wherein said modified 
tessellating is as follows: 

for all N x,y,z = v(x...x + 1, y...y + 1, z,..z + 1) 
1 0 compute octant o(pc, y, z) 

index polygon look-up table using o(x, y, z) 

add triangles corresponding to o(x, y, z\ offset by (pc, y, z). 

end 

123. An article of manufacture as in Claim 121 wherein said filtering comprises 

1 5 signal processing as follows: 

for fc= 1 : n 

for all e V {for each vertex} 
SXx,yyZ)4r-Q 

for all Vj e adjXVg) {for each vertex adjacent to Fi} 

20 S,(x,y,z)^S,+Vj 

end 

V,<^(l^a)V,^a(S,/\adj(y,)\) 
end 

end 

25 

124. A system for generating a smoothed surface representation of a three- 
dimensional voxel representation of a three-dimensional image comprising: 

a processor configured to: 

a. segment a three-dimensional image to provide a segmented 
30 voxelated three-dimensional signal; 

b. tessellate said segmented 3D image to provide a 3D triangulated 
surface representation; and 

c. filter said 3D representation resulting firom (b) to smooth said 3D 
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triangulated surface representation. 



125. A system as in Claim 124 wherein said processor tessellates as follows: 

for all N x,y,z = v(x...x + 1, y,.,y + 1, z.,.z + 1) 
5 compute octant y, z) 

index polygon look-up table using o(x, y, z) 

add triangles corresponding to o(x, y, z\ offset by (jc, y, z). 

end 

126. A system as in Claim 124 wherein said processor filters said 3D 
1 0 representation by signal processing as follows: 

for k—l:n 

for all K/ e V {for each vertex} 
Sfix,y,z)<r'Q 

for all e adj(y^) {for each vertex adjacent to Vi) 
15 S,ix,y,z)^S,+Vj 

end 

V, ^il^a)V,^aiS,l\adj(y,)\) 

end 
end 

20 

127. A method of estimating volumetric doubling time (DT) of an object 
which changes size, said object found in a maromalian host, comprising: 

a. obtaining a measurement of a change of a volxmietric characteristic 
of an object in a mammal over a time period, said object characterized as 

25 changing in size, over a time period, 

b. relating said change in said volumetric characteristic foimd in step 
(a) to said time period via an exponential size change model; and 

c. determining an estimated doubling time (DT) by comparing 
volume change over said time period. 

30 

128. A method as in Claim 122 wherein said measurement is obtained by two 
three-dimensional volumetric assessments of an existing object at a first time (fi) and 
again at later second time (ti). 
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129. A method as in Claim 128 wherein said exponential size change model is 

* 

wherein Vi and V2 are first and second volumetric assessments, respectively, and A is an 
exponential coefficient relating to the doubling time and 

DT 

130. A method as in Claim 129 wherein said doubling time (DT) is determined 
as follows: 

]n2-At 

131. A method as in Claim 127 which ftirther comprises combining a 
10 retardation factor (R) with said standard exponential size change model. 



15 



25 



132. A method as in Claim 127 wherein said measurements are obtained from 
three-dimensional images provided j&om signals bearing three-dimensional data produced 
firom scanning said host. 



133. A method as in Claim 132 wherein said image has been supersampled. 



134. A method as in Claim 132 wh^ein said image has been segmented. 



20 135. A method as in Claim 132 wherein said image has been supersampled and 

segmented. 



136. A method as in Claim 132 wherein said host is human and said object is a 

growth. 



137. A method as in Claim 127 and wherein said growth is a pulmonary 

nodule. 
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138. A method as in Claim 137 wherein said puhnonary nodxile is not greater 
than 10 mm in diameter. 

5 139. A method as in Claim 137 wherein said measurement of volume change is 

obtained from an object which is not detected at a first time (fi), but which is detected at a 
second time(/2). 

140- A method as in Claim 139 which comprises: 
10 (i) defining a minimum-detectable size dimension P ; and 

(ii) calculating doubling time using volume measurement a t2 
and a volume measurement derived from said minimum-detectable 
growth dimension jS, and the time period (A/^) between and t2. 

15 141, A method as in Claim 1 40 which comprises a determination of doubling 

time estimation as follows: 

142. A method as in Claim 127 wherein said volumetric characteristic is mean 
density and said change is change in mean density. 

20 

143. An article of manufacture for estimating volimietric doubling time (DT) of 
an object which changes size, said object found in a marmnalian host, comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: 
2S a. obtaining a measurement of a change of a volumetric characteristic 

of an object in a mammal over a time period, said object characterized as 
changing in size, over a time period, 

b. relating said change in said volumetric characteristic found in step 
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(a) to said time period via an exponential size change model; and 

c. detemiining an estimated doubling time (DT) by comparing 
volume change over said tune period. 



144. An article of manu&cture as in Claim 143 wherein said measurement is 
obtained by two three-dimensional volumetric assessments of an existing object at a first 
time (^i) and again at later second time (^2)- 



145. An article of manufacture as in Claim 144 whereiu said exponential size 
change model is 

wherein Vy and V2 are first and second volumetric assessments^ respectively, and A is an 
exponential coefficient relating to the doubling time and 

DT 

146. An article of manufacture as in Claim 145 wherein said doubling time 
(DT) is determined as follows: 

j^ ^j^ In 

■"ln(K,/F,) 

147. An article of manufacture as in Claim 143 which fiirther comprises 
combining a retardation factor (R) with said standard exponential size change model. 



148, An article of manufacture as in Claim 143 wherein said measurements are 
obtained firom three-dimensional images provided firom signals bearing three-dimensional 
data produced from scanning said host. 



149. An article of manufacture as in Claim 148 wherein said image has been 
supersampled. 
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150. An article of manufacture as in Claim 148 wherein said image has been 
segmented. 

151. An article of manufacture as in Claim 1 48 wherein said image has been 
5 supersampled and segmented. 

152. An article of manufacture as in Claim 143 wherein said host is human and 
said object is a groAvth. 

10 153. An article of manufacture as in Claim 143 and wherein said growth is a 

pulmonary nodiile. 

154. An article of man\rfacture as in Claim 153 wherem said pxilmonary nodule 
is not greater than 10 nun in diameter. 

15 

155. An article of manufacture as in Claim 143 wherein said measurement of 
volume change is obtained from an object which is not detected at a first time (^i), but 
which is detected at a second time(/2)- 

20 1 56. An article of manufacture as in Claim 1 55 which comprises: 

(i) defbiing a minimum-detectable size dimension j3 ; and 

(ii) calculating doubling time using volimie measurement a t2 
and a volume measurement derived from said minimum-detectable 
growth dimension J3, and the time period (A/^) between and ^2- 

25 

1 57. An article of manufacture as in Claim 1 56 which comprises a 
determination of doubling time estimation as follows: 
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In2-Ar, 



ln(F,/F^) 



158. An article of manufacture as in Claim 143 wherein said volumetric 
characteristic is mean density and said change is change in mean density. 

159. A system for estimating volumetric doubling time (DT) of an object which 
changes size, said object foxmd in a mammalian host, comprising: 

a processor configured to: 

a. obtain a measurement of a change of a volumetric characteristic of 
an object in a mammal over a time period, said object characterized as changing 
in size, over a time period, 

b. relate said change in said volumetric characteristic found in step 
(a) to said time period via an exponential size change model; and 

c. determine an estimated doubling time (DT) by comparing volume 
change over said time period. 

160. A system as in Claim 1 59 wherein said processor is further configured to 
obtain said measurement by two three-dimensional volumetric assessments of an existing 
object at a first time (A) and again at later second time (Ji). 

161. A system as in Claim 160 wherein said exponential size change model is 

wherein V\ and V2 are first and second volumetric assessments, respectively, and ;i is an 
exponential coefficient relating to the doubling time and 

DT 

162. A system as in Claim 161 wherein said doubling time (DT) is determined 
as follows: 
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In2-A/ 



163. A system as in Claim 159 wherein said processor is further configured to 
combine a retardation factor (R) with said standard exponential size change model. 



164. A system as in Claim 159 wherein said processor is further configured to 
obtain said measurement firom three-dimensional images provided fiom signals bearing 
three-dimensional data produced fiom scanning said host. 



165. A system as in Claim 164 wherein said image has been supersampled. 



166. A system as in Claim 164 wherein said image has been segmented. 



1 67. A system as in Claim 164 wherein said image has been supersampled and 
segmented. 



168. A system as in Claim 159 wherein said host is hxmian and said object is a 

growth. 



1 69. A system as in Claim 159 wherein said growth is a pulmonary nodule. 



1 70. A system as in Claim 169 wherein said pulmonary nodule is not greater 
than 1 0 mm hi diameter. 



171. A system as in Claim 159 wherein said measurement of volume change is 
obtained fiom an object which is not detected at a first time (/i), but which is detected at a 
second time(f2). 



172. A system as in Claim 171 wherein said processor is further configured to: 
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(i) define a mmimum-detectable size dimension j3 ; and 

(ii) calcxilate doubling time using volume measurement a tz and 
a volume measurement derived from said minimum-detectable 
growth dimrasion jS, and the time period (A/^) between /i and /2- 



173. A system £is in Claim 172 wherein said processor calculates a 
determination of doubling time estimation as follows: 



ln(K,/Fp) 

« 

10 1 74. A system as in Claim 1 59 wherein said volumetric characteristic is mean 

density and said change is change in mean density. 

175. A method for determining the time between scaiming examinations 
(InterScan interval) of an object which changes in size and is found in a host, comprising: 

15 a. developing a maximum magnitude of error (e for maximum 

percent error) between volumetric measurements derived from scanning and 
actual object size, said relationship specific to the nature of said object; 

b. determining a minimum reliably-detectable percent volume change 
(a) based on said maximum magnitude of ettor; 
20 c. establishing a relationship between magnitude of error and object 

size; 

d. deriving an expression using an exponential growth model in terms 
of change in volume and time needed to observe a change in volume; and 

e. based on the expression resulting from determining a time 
25 interval to observe a reliably-detectable percent change. 



176 
follows: 



A method as in Claim 175 whrarein said interscan interval is calculated as 
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^DT^ .ln(l + (a/100))/ln2 
wherein DTb is an estimated doubling time, and a is 2 - 6. 

177. A method as in Claim 175 wherein said object is a potentially cancerous 
growth and said InterScan interval is sufficiently short in duration to increase probability 
of early detection as cancerous. 

178. A method as in Claim 175 wherein said growth is a pulmonary nodule. 

179. A method as in Claim 178 wherein said nodule is not greater than 10 mm 
in diameter. 

1 80. A method for determining the time between scaiming examinations 
(InterScan interval) of an object which changes in size and is found in a host, comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: 

a. developing a maximum magnitude of error (€ for maximum 

percent error) between volimietric measurements derived from scanning and 
actual object size, said relationship specific to the nature of said object; 

b. determining a minimum reliably-detectable percent volume change 
(a) based on said maximum magnitude of error; 

c. establishing a relationship between magnitude of error and object 

size; 

d. deriving an expression using an exponential growth model in terms 
of change in volume and time needed to observe a change in volume; and 

e. based on the expression resulting from J, determining a time 
interval to observe a reliably-detectable percent change. 

181. An article of manufacture as in Claim 1 80 wherein said InterScan interval 
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is calculated as follows: 



r,- = Dr^ * ln(l + (a /lOO)) /ln2 



wherein DTb is an estimated doubling time, and a is 2 • € . 



5 



1 82. An article of manufacture as in Claim 1 80 wherein said object is a 



potentially cancerous growth and said InterScan interval is sufficiently short in duration to 
increase probability of early detection as cancerovis, 

183. An article of manufacture as in Claim 1 80 wherein said growth is a 
1 0 piilmonary nodule. 

1 84. An article of manufacture as in Claim 1 83 wherein said nodule is not 
greater than 10 mm in diameter. 

15 1 85. A system for determining the time between scanning examinations 

(InterScan interval) of an object which changes in size and is found in a host, comprising: 



a process configured to: 

a. develop a maximum magnitude of error (e for maximum percent 



error) between volumetric measurements derived from scaiming and actual object 



20 



size, said relationship specific to the nature of said object; 

b. determine a minimum reliably-detectable percent volume change 



(a) based on said maximimi magnitude of error; 

c, establish a relationship between magnitude of error and object size; 

d. derive an expression using an exponential grovsrtih model in terms 



25 



of change in volume and time needed to observe a change in volume; and 

e. determine a time interval to observe a reliably-detectable percent 



change based on the expression. 



186. A system as in Claim 185 wherein said iaterscan interval is calculated as 
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follows: 

t, ^DTj, •ln(H-(a/100))/ln2 
wherein DTd is an estimated doubling time, and a is 2 - £. 

187. A system as in Claim 185 wherein said object is a potentially cancerous 
growth and said InterScan interval is sufficiently short in duration to increase probabiUty 
of early detection as cancerous. 

188. A system as in Claim 185 wherein said growth is a pulmonary nodule. 

1 89. A system as in Claim 1 88 wherein said nodule is not greater than 10 mm 
in diameter. 

190. A method of determining volumetric change of size over a period of time 
comprising: 

a. determining a doubling time (DT) for a change in volumetric 
dimension of an object foimd in a host; 

b. establishing a relationship based on an exponential change model 
and said period of time between volumetric doubling time and the relative change* 
in volume; and 

c. determining said volumetric change over said period of time. 

191. A method as in Claim 190 wherein said period is one year (365.25 days) 
and said relationship is: 

^r^^ / hi2 • 365.25 days x ^ 
VGI- exp ( —) -1 

1 92. An article of maniifacture for determining volumetric change of size over 
a period of time comprising: 
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a machine readable-medivim containing one or more programs which 
when executed implement the steps of; 

a* deteraiining a doubling time (DT) for a change in volumetric 
dimension of an object found in a host; 

b. establishing a relationship based on an exponential change model 
and said period of time between volumetric doubling time and the relative change 
in volvime; and 

c. determining said volumetric change over said period of time. 

193. An article of manufacture as in Claim 192 wherein said period is one year 
(365 ;25 days) and said relationship is: 

^ hi2 -36 5,25 days , 
FG/= exp ( —) -1 



194. A system for determining volumetric change of size over a period of time 
comprising: 

a processor configured to: 

a. determine a doubling time (DT) for a change in volumetric 
dimension of an object found in a host; 

b. establish a relationship based on an exponential change model and 
said period of time between volumetric doubling time and the relative change in 
volume; and 

c. determine said volumetric change over said period of time. 



195. A system as in Claim 194 wherein said period is one year (365.25 days) 
and said relationship is: 

1/^7- / In 2 > 365.25 days ., . 
VGI= exp ( ) -1 

1 96. A method of analyzing three-dimensional scan-produced images.of a 
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region of interest (ROI) taken at least two separate times, comprising three-dimensionally 
registering two 3D ROI images taken at separate times. 

197. A method as in Claim 196 wherein said three-dimensional registration 
5 comprises locating the center of mass (COM) in each said image of said ROI, and 
aligning same for said analysis. 



198. A method as in Claim 197 wherein said locating comprises iteratively 
determining center of mass of an iteratively refined ROI. 

10 

199. A method as in Claim 198 wherein a current COM, Cw, is used to replace 
a previous COM, Cc, and said ROI is recomputed until the distance between the COM is 
not greater than some pre-selected smaH distance. 



15 200. A method as in Claim 199 wherein said iterative COM determination 

comprises: 

Select an initial location for the COM: Q = v(Xi,yi,Zi) 
Select the desired size of the ROI: S = x,y,z 
AC = oo 

20 Cc = Ci 

whUe (AC > e) 

Compute new ROI bounds based on Cc and S 

C,k(x) = fwioo/'wooo 

Cm(y) = moio/mooo 
25 Cm(z) = wooi/'Wooo 

AC = D(Cm. Co) 
Cc Cm 

end 

where DiVa^Vt) is the Euclidean distance function between voxels Va and Vi,. 

30 

20 1 . A method as in Claim 1 96 wherein said three-dimensional registering 
comprises 3D correlation. 
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202. A method as in Claim 201 wherein said correlation comprises: 

a. selecting two ROI images; 

b. specifying the limits of a search region; 

c. translating a second ROI to all locations in the search region; and 
5 d. computing the 3D correlation between the two ROIs to determine 

the greater value of a match metric. 

203, A method as in Claim 202 which is implemented by: 

Select two ROIs: ROIa, ROh 
1 0 Select the extent of the search region in each dimension: S = {x, z) 

for A: = -iSr : iSz {for each translation in search region} 
for j = 'Sy:Sy 
for i — -Sx:Sx 

15 c(i, j\ k) = ROI^ ^ ROJb 

if eft/ k)>M 

M=c(i,j\ k) 

L'^ihj, k) 
end 

20 end 

end 
end 

wherein T(ROI, j, k,) specifies the translation of an ROI coordinate system by x=7, 
25 y—J,z=k^ and the star(^) operator is a three-dimensional correlation as follows: 

where in each iteration the correlation of the two ROIs at the current translation of ROI is 
computed until the value of the correlation, is greater than the previous m aximum 
value My so that the value and the location of the better match , L, are recorded. 



30 



204. A method as in Claim 196 which further comprises three-dimensional 
weighting said ROI to reduce unwanted confounding structures found in said ROL 
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205. A method as in Claim 204 wherein said three-dimensional weighting is 
implemented using one of: (i) a Gaussian window, (ii) a triangular radial window, and 
(iii) a rectangular radial window. 

5 206. A method as in Claim 204 which further comprises imposing the method 

of moment analysis to provide density distribution of said weighted ROI. 

207. A method as in Claim 204 which further comprises estimating doubling 
time (DT) of a growth in said ROI based on change in mean density of said ROI. 

10 

208. A method as in Claim 207 wherein the growth image is Gaussian- 
weighted and said estimation comprises an iterative technique to select an appropriate 
standard deviation, a, for weighting. 



15 209. A method as in Claim 208 wherein said iterative technique is as follows: 

Select an ROI, ROI, and a standardized mean density value, mds, 

5 = 0.1 

mdc-0 

20 while done 9^ 1 

Generate Gaussian Weighting function Wc using ac 
WROIcQc, y, z)^w^x, y, z) • ROJQc, y, z) 

md,^(x^' E^^o Z^i WROI,ix,y,z))/iM'N^L) 

8 = {mdj— mdc) 
25 ifl(5|>0 

Ob = Oc 

if 1^1 < 6: 

done = 1 
else 

30 (Tc^CTc'^^ 

end 

else 

if |i5|<£: 
done = 1 

35 else 

s = s/2 
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end 
end 

end 



for a given standardized mean density, mdsy and stopping criteria, € . 

210. A method as in Claim 209 wherein doubling time estimation comprises: 

Select two ROIs, ROIj, ROh, at t} and t2, respectively, 
1 0 Align {ROIu ROh) using correlation-based region matching. 

Determine a based on ROIi using iterative sigma determination. 
Generate Gaussian weighting function w using a 
WROhipc, z) = w{x, z) • ROIi(x, y, z) 
WROhix, y, z) = >v(x, y, z) • ROhix, y, z) 

15 md, =(Z^' ZiL'o Zfco WROI,(x,y,z))KM^N'L) 

md2^(X^' E^l WROI^{x,y,z))l(M'N^L) 

In2-A/ 

DTj. = — 

hiQndj/md^) 

211, An article of manufacture for analyzing three-dimensional scan-produced 
20 images of a region of interest (ROI) taken at least two sepamte times, comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: three-dimensionally registering two 3D 
ROI images taken at separate times. 

25 2 1 2. An article of manufacture as in Claim 211 wherein said three-dimensional 

registration comprises locating the center of mass (COM) in each said image of said ROI, 
and aligning same for said analysis. 

213. An article of manufacture as in Claim 212 wherein said locating comprises 
30 iteratively detennining center of mass of an iteratively refined ROI. 
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214. An article of manufacture as in Claim 213 wherein a cvirrent COM, Cm, is 
used to replace a previous COM, Cc, and said ROI is recomputed until the distance 
between the COM is not greater than some pre-selected small distance. 

5 215. An article of manufacture as in Claim 214 wherein said iterative COM 

determination comprises: 

Select an initial location for the COM: C/ = v(xi,yi,Zi) 
Select the desired size of the ROI: S = x,y,z 
AC=oo 

10 Cc = C, 

while (AC >e) 

Compute new ROI bounds based on Cc and S 

Cm(x) = /Wioo/WoOO 

Cm(y) = moxo/mm 
1 5 Cm(z) = wooiZ/wooo 

AC^DiCm^Cc) 
Cc Cm 

end 

where DiVa^Vt) is the Euclidean distance function between voxels Va and Vb^ 

20 

216. An article of manufacture as in Claim 211 wherein said three-dimensional 
registering comprises 3D correlation. 

217. An article of manufacture as in Claim 216 wherein said correlation 
25 comprises: 

a. selecting two ROI images; 

b. specifying the limits of a search region; 

c. translating a second ROI to all locations in the search region; and 

d. computing the 3D correlation between the two ROIs to determine 
3 0 the greater value of a match metric. 

218. An article of manufacture as in Claim 216 which is implemented by: 
Select two ROIs: ROIa, ROh 

Select the extent of the search region in each dimension: S = {x, z) 
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M=0 

for k—'Sz: Sz {for each translation in search region} 
for J = -Sy: Sy 
for ! = 

5 c(hj\ k) = ROh ^ ROh 

ilc(i,lk)>M 

M^cihj, ky 

end 

10 end 

end 
end 

wherein T(ROI, h j> h) specifies the translation of an ROI coordinate system by x^i, 
1 5 yz=j^z=k^ and the star(T*r) operator is a three-dimensional correlation as follows: 

X y 2 

where in each iteration the correlation of the two ROIs at the current translation of ROI is 
computed vmtil the value of the correlation, c, is greater than the previous maximum 
value M, so that the value and the location of the better match, are recorded. 

20 

21 9. An article of manufacture as in Claim 211 which further comprises three- 
dimensional weighting said ROL 

220. An article of manufacture as in Claim 2 19 wherein said three-dimensional 
25 weighting is implemented using one of: (i) a Gaussian window, (ii) a triangular radial 

window, and (iii) a rectangular radial window. 

221 . An article of manufacture as in Claim 219 which further comprises 
imposing the method of moment analysis to provide density distribution of said weighted 

30 ROI. 



222. An article of manufacture as in Claim 219 which further comprises 
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estimatiiig doubling time (DT) of a growth in said ROI based on change in mean density 
ofsaidROL 



223. An article of manufacture as in Claim 222 wherein the growth image is 
5 Gaussian-weighted and said estimation comprises an iterative technique to select an 
appropriate standard deviation, a, for weighting. 



224. An article of manufacture as in Claim 223 wherein said iterative technique 
is as follows: 

1 0 Select an ROI, ROI^ and a standardized mean density value, mds^ 

mdc = 0 
while done ^ 1 

1 5 Generate Gaussian Weighting function Wc using Oc 

WROIcQc, z) = Wc(x:, z) • ROI{x, y, z) 

md, = (S^^ S^to Z^ri WROI, (X, y, z)) /(M^N^ L) 

S^imds—mdc) 
if |<51>0 

20 ab = (Tc 

if\d\<e 

done = 1 
else 

25 end 

else 

if |<5|<6: 
done = 1 

else 

30 5 = s/2 

end 
end 



35 



end 



for a given standardized mean density, mdg^ and stopping criteria, €. 
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225. An article of manufacture as in Claim 224 wherein doubling time 



estimation comprises: 



10 



5 





226. A system for analyzing three-dimensional scan-produced images of a 
1 5 region of interest (ROI) taken at least two separate times, comprising: 

a processor configured to three-dimensionally register two 3D ROI images 
taken at separate times. 

227. A system as in Claim 226 wherein said processor is further configured to 
20 locate the center of mass (COM) in each said image of said ROI, and align same for said 

analysis. 

228. A system as in Claim 227 wherein said processor is configured to locate 
the center of mass by determining the center of mass (COM) of an interatively refined 



229, A system as in Claim 228 wherein a current COM, Ca,, is used to replace a 
previous COM, Cc, and said ROI is recomputed until the distance between the COM is 
not greater than some pre-selected small distance. 



25 



ROL 



30 



230. A system as in Claim 229 wherein said iterative COM determination 
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comprises: 

Select an initigd location for the COM: C/ = y(Xi,yi,Zi) 
Select the desired size of the ROI: S = x,y,z 
AC-oo 

Cc = Ci 
while (AO e) 

Compute new ROI bounds based on Cc and S 

Cm(x) = mioo/moQo 

Cjn(y) = mio/mm 
Cm(z) = moQi/mQQQ 

AC = i)(C^, Cc) 

Cc C/n 

end 

where D(Va>Vb) is the Euclidean distance function between voxels Va and Vb- 



231, A system as in Claim 226 wherein said processor is configured to register 
by a 3D correlation. 



232. A system as in Claim 23 1 wherein said correlation comprises: 

a. selecting two ROI images; 

b. specifying the limits of a search region; 

c. translating a second ROI to all locations in the search region; and 

d. computing the 3D correlation between the two ROIs to determine 
the greater value of a match metric. 



233. A system as in Claim 23 1 which is implemented by: 
Select two ROIs: ROIa, ROh 

Select the extent of the search region in each dimension: S = {x^ z) 
for k-'Sz'. Sz {for each translation in search region} 

for J—'SylSy 

for i--Sx:Sx 

c{iJ,k)^ROh^ ROh 
ifc{ij\k)>M 

M^c(J,j\ k) 

L^ihj, k) 

end 
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end 

end 

end 

5 wherein T(ROI, ij, k,) specifies the translation of an ROI coordinate system by x=z, y^j, 
z=^, and the star(T*r) operator is a three-dimensional correlation as follows: 

X y z 

wherein each iteration the correlation of the two ROIs at the cmrent translation of ROI is 
computed until the value of the correlation, c, is greater than the previous maximum 
10 value M, so that the value and the location of the better match, L, are recorded. 

234. A system as in Claim 226 wherein said processor is further configured to 
three-dimensionally weigh said ROI. 

15 235. A system as in Claim 234 wherein said processor is configured to weigh 

said ROI by using one of: (i) a Gaussian window, (ii) a triangular radial window, and (iii) 
a rectangular radial window. 

236. A system as in Claim 234 wherein said processor is configured to impose 
20 the method of moment analysis to provide density distribution of said weighted ROI. 

237. A system as in Claim 234 wherein said processor is configured to estimate 
doubling time (DT) of a growth in said ROI based on change in mean density of said 
ROI. 

25 

238. A system as in Claim 237 wherein the growth image is Gaussian- 
weighted and said estimation comprises an iterative technique to select an appropriate 
standard deviation, c, for weighting. 
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239. A system as in Claim 238 wherein said iterative technique is as follows: 
Select an ROI, ROI, and a standardized mean density value, mds^ 
s = QA 

5 mdc = 0 

while done 9^ 1 

Generate Gaussian Weighting function Wc using Gc 
WROIc(x, z) = wJix, y, z) - ROIipc, y, z) 

md.^iX'^' 2^ WROI,{x,y,z))l{M^N^D 

10 d-(rnds-mdc) 

if|(51>0 

if\d\<€ 
done = 1 
IS else 

end 

else 

il\d\<e 

20 done — 1 

else 

s = sl2 

end 

25 end 

end 

for a given standardized mean density, mds, and stopping criteria, € . 



240. A system as in Claim 239 wherein doubling time estimation comprises: 

30 Select two ROIs, ROIj, ROI29 at t] and /2> respectively. 

Align (ROIjy ROI2) using correlation-based region matching. 
Determine <r based on ROIi using iterative sigma determination. 
Generate Gaussian weighting function w using a 
WROIj(x, y, z) = w(x, y, z) • ROIi(x, y, z) 

35 WROhix, y, z) = y, z) • ROhix, y, z) 

md,^(X^' E^o Zz^o WROI,(x,y,z))/(M^N'L) 
mrf2=(z^o' Z^LI WROI^ix,y,z))l{M N'L) 
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In2-Ar 



24 1 . A method of estimating curvature of an object which has been scanned to 
provide a signal of a three-dimensional representation of said object comprising: 

(a) providing a three-dimensional triangularly tessellated 
representation of said object; 

(b) determining the surface normal to all triangles; 

(c) based on said surface noraials result from step (b), calculating 
vertex surface normal at each vertex; 

(d) determining the angular difference between the surface normal at 
each vertex (designated home vertex, and the vertex normals of all vertices 
adjacent to said home vertex; and 

(e) estimating curvature at each vertex of said object based on angular 
differences resulting from step (d). 



242. A method as in Claim 241 which further comprises detecting triangles 
result, from step (a) 



243. A method as in Claim 242 wherein said triangle detecting is as follows: 
VKi G V 

V(F^, Vb G adj (Vi) {For each pair of vertices adjacent to V/} 
if Fb e adj (Va) {Ifadj (Va) contains Vt} 

end 

244, A method as in Claim 241 wherein said angular difference of step (d) is 
determined as follows: 



(l>^ = cos 



245. A method as in Claim 241 wherein said ctiirvature estimate includes 
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principal curvatures, ku and determined as follows: 

for all Vg € K(for each vertex) 
for allF^ea^^'CF;) 

0,=cos-'((/),-0J 

end 

Ki(7J=max(0) 
fC2(F,)=min(0) 
end 

246. A method as in Claim 241 wherein said curvature estimate is a Gaussian 
curvature wherein 

K =K^K2 • 

247. A method as in Claim 241 wherein said curvature estimate as a mean 
curvature, Hy wherein 

jTZ = . 

2 

248. A method as in Claim 241 wherein said curvature estimate is an average 
normal curvature estimate as follows: 

for all G F (for each vertex) 

(compute curvature at vertex as average of angular differences) 
for all e adJiV.) 

0„=COS-H(-^,-<^a)/(|^,|K|)) 

end 
end 

249. A method as in Claim 241 wherein all steps (b) - (e) are implemented as 
follows: 



180 



wo 01/78005 



PCTAJSOl/11820 



for =V^ (for each vertex) 
for (V^,n)sadjXV,) 
ifV,^adJ{V,) 
z=>3Tj={V,,V^,V^}iV,is2i memberof trianglcT;-) 

end 

for Tj =To :T^ (for each triangle Tj) 

XI/ J =(/a X ib X ib (calculate triangle surface normal) 
end 

®i =(SJ^o^^)'^(''^~^)(^^^^^^^*® vertex surface normal) 
end 

for V. =^Vq : V„ (for each vertex) 
for V,^adf(V,) 

0,=cos-'^((^,.0J/(|^,||^J)| 

= ^2^^^ J /i\^e^ (calculate curvature at vertex) 

end 
end 

for T.^To :7; (for each triangle) 

Cj,. ={Cy^ + ^vb)^^ (calculate curvature on triangle) 
end 

250. An article of manufacture for estimating curvature of an object which has 
been scanned to provide a signal of a three-dimensional representation of said object 
comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: 

(a) providing a three-dimensional triangularly tessellated 
representation of said object; 

(b) determining the surface normal to all triangles; 

(c) based on said surface normals result &om step (b), calculating a 
vertex surface normal at each vertex; 

(d) determining the angular difference between the surface normal at 
each vertex (designated home vertex, Fi) and the vertex normals of all vertices 
adjacent to said home vertex; and 

(e) estimating at each vertex curvature of said object based on angular 
differences resulting from step (d). 
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25 1 . An article of manufacture as in Claim 250 which further comprises 
detecting triangles result, from step (a) 



252, An article of manufacture as in Claim 25 1 wherein said triangle detecting 
is as follows: 

V(Kfl, e adj (Vi) {For each pair of vertices adjacent to V/} 
if Fb G adj (Va) {IfadJ (Va) contains Fi,} 

end 

253. An article of manufacture as in Claim 250 wherein said angular difference 
of step (d) is determined as follows: 



0^ = cos'^ 



{M<f>a\) 

254. An article of manufacture as ia Claim 250 wherein said curvature estimate 

includes principal curvatures, ku and K2, determined as follows: 

for all e F(for each vertex) 
for allF^ea^/'C'^) 

0„ =cos (0, 

end 

K:,(F,)=max(e) 
Kj(F,)=mm(0) 
end 



255. An article of manufacture as in Claim 250 wherein said curvature estimate 
is a Gaussian curvature AT, wherein 

256. An article of manufacture as in Claim 250 wherein said curvature estimate 
as a mean cxxrvature, if, wherein 
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2 

257. An article of manufacture as in Claim 250 wherein said curvature estimate 
is an average normal curvature estimate as follows: 



for all € F (for each vertex) 

(compute curvature at vertex as average of angular differences) 
for aU e adj{V,) 

0,=cos-H(0r0J/(l</'/|W)) 

end 

end 

258. An article of manufacture as in Claim 250 wherein all steps (b) - (e) are 
implemented as follows: 

for V^-V^ : V„ (for each vertex) 
for (V,.K)^adjXK) 
if V,eadj(VJ 

=>3 Z}={F;,^;,F^}(F,isa member of triangle r,) 

end 

for Tj =To:T^ (for each triangle Tj ) 

xffj =(za X ib )/ ia x ib (calcxilate triangle surface normal) 
end 

=(X>o^y)^('^~^)(^^^^^^*® vertex siuf ace normal) 

end 

for V^^Vq : V„ (for each vertex) 
for V,^adJ(V,) 

0,=cos-^((^,.(^.J/(0,||^J)j 

end 
end 

for T, =Tg :T„ (for each triangle) 

^Ti "^ypva '^^vb ^vbT ^ (calculate curvature on triangle) 
end 
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259. A system for estimating ciirvature of an object which has been scanned to 
provide a signal of a three-dimensional representation of said object comprising: 

a processor configured to: 

(a) provide a three-dimensional triangularly tessellated representation 
5 of said object; 

(b) determine the surface normal to all triangles; 

(c) calculate a vertex surface normal at each vertex based on said 
surface normals; 

(d) determine the angular difference between the surface normal at 
10 each vertex (designated home vertex, Vi) and the vertex normals of all vertices 

adjacent to said home vertex; and 

(e) estimate curvature at each vertex of said object based on angular 
differences resulting from step (d). 

15 260. A system as in Claim 259 wherein said processor is further configured to 

detect triangles result based upon the three-dimensional triangularly tessellated 
representation of said object. 



261 . A system as in Claim 260 wherein said processor detects said triangles as 
20 follows: 

\fVi e V 

V(Fa. Vb G adj (Fi) {For each pair of vertices adjacent to Vj 

if Kb G adJ (Va) {if adj (Va) contains Vt} 

T^(TU{Vi,Va.Vt}) 
25 end 



262. A system as in Claim 259 wherein said processor is configured to 
determine said angular difference as follows: 



<lf„ = cos' 
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263. A system as in Claim 259 wherein said cmrvature estimate includes 
principal cmrvatures, and /ca, detemiined as follows: 

for all Vf € K(for each vertex) 
for siHV^Gadj(y,) 



end 

Kj(F'^)=max(0) 

K2(F;)=min(0) 
end 



264. A system as in Claim 259 wherein said processor is conJBgured to estimate 
said curvature as a Gaussian curvature wherein 

265. A system as in Claim 259 wherein said processor is configured to estimate 
said curvature as a mean curvature, H, wherein 

2 

266. A system as in Claim 259 wherein said processor is configvired to estimate 
said curvature as an average normal curvature as follows: 

for all e F (for each vertex) 

(compute curvature at vertex as average of angular differences) 

for all F; g adjXV.) 

e,=cos-K(<|»r0J/(|0/^»|)) 

end 
end 

267. A system as in Claim 259 wherein said processor is configured to 
implement all steps (b) — (e) as follows: 
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for :V„ (for each vertex) 

for (V,,V,)eadjXV,) 

=>3 Tj ={V,,V„,V^}(V, isa memberof triangleT}) 

end 

for Tj ^Tq (for each triangle Tj) 

y^j =(ia X ib )/j/a x ib (calculate triangle surface normal) 
end 

Oj =(^J^VAy)/(m=l)(calculate vertex surface normal) 
end 

for V^—Vq :V„ (for each vertex) 
for V^^adjXV,) 

0,=cos-^((^,-(^J/(|<^,||i^J)| 

Cyi =^y^!^ Qg^l ^Brt (calculate cxirvature at vertex) 

end 

end 

for 7;. =7^0 (for each triangle) 
Ti={ya,V„VJ 

^Ti "^^Pva ^^vb '^^vb)^^ (calculate curvature on triangle) 
end 



268. A method of providing a curvature metric analysis of a three-dimensional 
object in a host which has been scanned to provide signed corresponding to said objects 
comprising: 

(a) providing a measurable three-dimensional representation of an 
object, said representation including a surface; 

(b) estimating a surface curvature at a number of discrete points on 
said surface to provide a number of surface curvature estimates; and 

(c) determining a frequency distribution of curvature estimates 
resulting from said surface. 



269. A method as in Claim 268 which further comprises determining a variance 
of said distribution. 

270. A method as in Claim 268 which ftirther comprises deterniining a range of 
said distribution. 

186 



wo 01/78005 



PCTAJSOl/11820 



10 



271. A method as in Claim 268 which further comprises determining a mean of 
said distribution. 



5 272. A method as in Claim 268 which further comprises determining a 

minimum of said distribution. 



273. A method as in Claim 268 which fixrther comprises determining a 
maximum of said distribution. 



274. A method as in Claim 268 wherein said object is a growth in a human. 



275. A method as in Claim 274 wherein said growth as a puhnonary nodule. 



1 5 276. A method as in Claim 268 wherein said nodule is not greater than 1 0 mm 

in diameter. 



277. An article of manufacture for providing a curvature metric analysis of a 
three-dimensional object in a host which has been scanned to provide signal 
20 corresponding to said objects comprising: 

a machine readable-medium containing one or more programs which 
when executed implement the steps of: 

(a) providing a measurable three-dimensional representation of an 
object^ said representation including a surface; 
25 (b) estimating a surface curvature at a number of discrete points on 

said surface to provide a number of surface curvature estimates; and 

(c) determining a frequency distribution of curvature estimates 
resulting from said surfsice. 



30 



278. An article of manufacture as in Claim 277 which further comprises 

187 



wo 01/78005 



PCTAJSOl/11820 



detennining a variance of said distribution. 

279. An article of manufacture as in Claim 277 which further comprises 
determining a range of said distribution. 

5 

280. An article of manufacture as in Claim 277 which further comprises 
determining a mean of said distribution. 

28 1 . An article of manufacture as in Claim 277 which further comprises 
1 0 determining a minimum of said distribution. 

282. An article of manufacture as in Claim 277 which further comprises 
determining a maximum of said distribution. 

1 5 283 . An article of manufacture as in Claim 277 wherein said object is a growth 

in a human. 

284. An article of manufacture as in Claim 283 wherein said growth as a 
pxilmonary nodule. 

20 

285. An article of manufacture as in Claim 277 wherein said nodxile is not 
greater than 10 mm in diameter. 

286. A system for providing a curvature metric analysis of a three-dimensional 
25 object in a host which has been scanned to provide signal corresponding to said objects 

comprising: 

a processor configured to: 

(a) provide a measmrable three-dimensional representation of an 
object, said representation including a surface; 
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(b) estimate a surface curvature at a number of discrete points on said 
surface to provide a number of surface curvature estimates; and 

(c) deteraiine a frequency distribution of curvature estimates resulting 
from said surface. 

287. A system as in Claim 286 v^herein said processor is further configured to 
determine a variance of said distribution. 

288. A system as in Claim 286 wherein said processor is further configured to 
determine a range of said distribution. 

289. A system as in Claim 286 wherein said processor is further configured to 
determine a mean of said distribution. 

290. A system as in Claim 286 wherein said processor is further configured to 
determine a minimimi of said distribution. 

291 . A system as in Claim 286 wherein said processor is further configured to 
determine a maximum of said distribution. 

292. A system as in Claim 286 wherein said object is a grovsrth in a human. 

293. A system as in Claim 292 wherein said growth as a pulmonary nodule. 

294. A system as in Claim 286 wherein said nodxile is not greater than 10 mm 
in diameter. 
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